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FOREWORD 


The  idea  for  this  report  originated  during  a  routine  review  of  several  water 
supply  modeling  studies.  These  studies  were  all  concerned  with  the  impacts  of  a 
proposed  increase  in  groundwater  withdrawals  from  the  San  Andrcs-Glorieta  aquifer 
in  northwestern  New  Mexico.  As  the  review  progressed  it  became  apparent  that  the 
drawdowns  predicted  by  the  various  studies  differed  dramatically,  even  though  they 
all  relied  on  the  same  data  base  and  computer  model.  This  prompted  a  number  of 
questions  -  Why  were  the  modeling  results  so  different?  Is  such  variabilty  typical? 
Can  any  conclusions  about  groundwater  modeling  in  general  be  drawn  from  this 
example? 

This  report  addresses  the  questions  which  grew  out  of  the  San  Andrcs-Glorieta 
modeling  review.  It  also  proposes  some  general  modeling  guidelines  which  are 
intended  to  deal  with  the  uncertainties  and  ambiguities  inherent  in  the  modeling 
process.  The  report  is  designed  for  use  in  training  courses  and  workshops  and 
so  is  written  in  a  tutorial  style.  Report  preparation  was  funded  under  Contract 
No.  DACW05-83-P-1410  from  the  U.S.  Army  Corps  of  Engineers,  Hydrologic 
Engineering  Center,  Davis,  CA.  95616. 
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1.  INTRODUCTION 


Over  the  last  lew  decades  computerized  groundwater  models  have  moved  from 
the  research  laboratory  into  the  offices  of  consulting  hydrologists  and  government 
planners.  The  wide  accessibility  of  computer  models  and  the  equipment  needed 
to  run  them  has  brought  about  dramatic  changes  in  the  way  groundwater  studies 
are  conducted.  This  is  particularly  true  in  water  supply  investigations.  Computer- 
generated  maps  are  commonly  found  in  water  supply  reports  and  computer  results 
are  frequently  entered  as  evidence  in  court  proceedings  dealing  with  water  rights. 

Computer  modeling  has  undoubtedly  helped  to  make  water  supply  planning 
more  reliable.  Ihit  the  widespread  use  of  computerized  technology  has  been,  in  some 
ways,  a  mixed  blessing.  The  highly  quantified  and  voluminous  outputs  generated  in 
modeling  studies  tend  to  give  the  impression  that  a  problem  is  better  understood 
than  it  really  is.  Sometimes  computer  results  can  actually  obscure  the  key  fac¬ 
tors  influencing  groundwater  flow  patterns,  particularly  if  these  results  are  poorly 
presented.  Computer  modeling  is  a  powerful  tool  but  it  needs  to  be  used  with 
discretion.  In  particular,  the  uncertainties  associated  with  the  modeling  process 
should  be  recognized  and  honestly  acknowledged. 

Perhaps  the  best  way  to  examine  the  role  of  uncertainty  in  a  modeling  study  is 
to  review  the  technical  decisions  which  must  be  made  when  a  model  is  formulated 
and  applied.  Important  assumptions  are  required  at  every  stage  of  the  process — 
when  the  model’s  equations  are  derived,  when  a  solution  procedure  is  selected, 
and  when  inputs  are  estimated.  Many  of  these  assumptions  are  ultimately  based 
on  subjective  interpretation  of  limited  amounts  of  ambiguous  held  data.  Different 
interpretations  lead,  of  course,  to  different  predictions,  making  the  modeling  process 
more  dependent  on  the  judgement  of  the  individual  modeler  than  is  generally 
recognized. 

This  report  is  concerned  primarily  with  the  technical  choices  which  influence 
the  predictions  made  by  groundwater  models.  In  order  to  keep  the  discussion 
specific,  attention  is  focused  on  a  typical  modeling  case  study.  The  particular 
example  selected  involves  an  application  by  Plains  Electric  Generation  and  Trans¬ 
mission  Cooperative,  Inc.  (Plains  Electric),  a  public  utility,  for  permission  to  pump 
groundwater  from  the  San  Andres-Glorieta  aquifer  in  northwestern  New  Mexico. 
As  part  of  its  application,  Plains  Electric  introduced  model  results  which  suggested 
that  the  impacts  of  its  proposed  pumpage  plan  would  be  acceptable.  Local  water 
users  responded  by  sponsoring  several  independent  modeling  studies.  These  studies 
all  used  the  same  groundwater  model  (developed  by  the  U.S.  Geological  Survey) 
but  *hc  inputs  and  general  approaches  to  model  application  differed  considerably. 
Some  of  the  modeling  studies  predicted  rapid  dewatering  in  the  vicinity  of  the 
Plains  well  field  while  others  predicted  only  minor  effects  on  local  water  levels. 
Overall,  the  introduction  of  computer  modeling  seemed  to  create  more  confusion 
than  it  dispelled. 


-  The  ease  study  analysis  presented  in  this  report  was  undertaken  with  the  belief 
that  the  San  Andres-Clorieta  experience  might  reveal  some  important  insights  about 
the  basic  strengths  and  weaknesses  of  modeling  technology.  Since  the  report  is 
designed  to  be  used  in  groundwater  training  courses,  it  is  written  in  a  tutorial  style 
which  focuses  on  general  principles  rather  than  technical  details.  The  references 
cited  throughout  the  text  provide  background  information  on  particular  subject 
areas  which  some  readers  may  wish  to  cover  in  more  depth. 

The  next  chapter  of  the  report  presents  background  material  on  the  San 
Andres-Clorieta  aquifer  system  and  briefly  summarizes  the  management  problem 
considered  in  the  case  study.  Chapter  3  reviews  basic  groundwater  modeling 
concepts,  with  emphasis  placed  on  the  technical  decisions  which  must  be  made 
in  a  typical  groundwater  impact  study.  Chapter  4  describes  the  various  modeling 
approaches  used  in  the  case  study  and  considers  how  each  modeler  dealt  with  the 
technical  issues  identified  in  Chapter  3.  The  report  concludes  with  a  suggested  set 
of  modeling  guidelines  for  groundwater  impact  investigations.  These  guidelines, 
which  are  based  both  on  the  case  study  and  on  general  principles,  attempt  to 
address  the  issues  of  data  interpretation  and  model  accuracy  in  a  systematic  way. 


2.  WATER  SUPPLY  IN  THE  SAN  ANDRES-GLORIETA  AQUIFER 


2.1  THE  MANAGEMENT  PROBLEM 

The  San  Andres-Glorieta  aquifer  is  part  of  the  San  Juan  basin,  a  layered 
geological  system  which  extends  throughout  much  of  northwestern  New  Mexico. 
The  portion  of  this  system  of  most  interest  here  lies  between  Gallup  and  Grants 
along  the  northern  flank  of  the  Zuni  mountains  (sec  Figure  2-l).  The  San  Andres- 
Glorieta  is  in  many  ways  a  typical  aquifer  of  the  Colorado  plateau.  It  is  a  sandstone- 
limestone  formation,  confined  by  adjoining  beds  of  shale  and  clay,  which  is  recharged 
from  an  outcrop  area  in  the  Zuni  range.  The  aquifer  was  once  the  source  of  many 
springs  and  flowing  wells  but  water  levels  have  gradually  dropped  in  recent  years 
and  most  groundwater  must  now  be  pumped. 

Groundwater  development  in  New  Mexico  is  closely  regulated  by  the  state, 
particularly  in  areas  which  arc  declared  to  be  “underground  water  basins”  by  the 
State  Engineer.  The  study  area  shown  in  Figure  2-1  covers  two  declared  basins, 
the  Gallup  basin  to  the  northwest  and  the  Bluewater  basin  to  the  southeast,  which 
are  hydrologically  continuous  but  are  separated  for  jurisdictional  purposes  by  the 
Continental  Divide.  In  1980  Plains  Electric  filed  an  application  to  withdraw  water 
from  a  well  field  located  near  the  existing  Shell  refinery  site  at  Ciniza.  This  water 
is  to  be  pumped  exclusively  from  the-San  Andres-Glorieta  formation.  Local  water 
users  who  believe  the  Plains  proposal  will  have  adverse  impacts  on  their  water 
supply  have  filed  protests  with  the  State  Engineer.  The  pumping  schedule  proposed 
by  Plains  Electric  is  keyed  to  projected  power  plant  operations  and  averages  about 
3600  acre-feet/year  over  the  38  year  plant  life  (see  Figure  2-2).  Since  this  represents 
an  increase  in  total  annual  pumpage  of  approximately  150  percent  above  1982  levels 
it  is  not  surprising  that  local  residents  are  concerned. 

The  water  supply  management  problem  posed  by  the  Plains  Electric  applica¬ 
tion  ultimately  involves  a  decision  to  accept,  reject,  or  add  stipulations  to  the  Plains’ 
proposal.  This  decision  will  be  based  in  part  on  an  evaluation  of  the  hydrologic 
effects  of  additional  pumpage.  If  projected  water  level  declines  are  moderate, 
remedial  measures  such  as  the  lowering  of  existing  well  pumps  may  be  sufficient  to 
protect  existing  water  rights.  If  projected  declines  are  more  severe,  such  measures 
may  not  be  adequate  and  the  pumping  plan  will  have  to  be  modified. 
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FIGURE  2-2 

PROJECTED  WATER  DEMANDS  FOR  THE  PLAINS  ELECTRIC  WELL  FIELD 


2.2  GEOIIYDROLOGY  OF  THE  SAN  ANDRES-GLORIETA  AQUIFER 


The  geohydrology  of  the  San  Andres-Gloricta  aquifer  has  been  briefly  reviewed 
in  most  of  the  modeling  studies  of  the  region  (see,  for  example,  Geohydrology, 
1982,  Gcotrans,  1982,  and  HEC,  1982).  Some  of  the  more  important  primary 
references  include  Shomaker  (1971),  Baars  and  Stephenson  (1977),  and  Read  and 
Wanck  (1961).  The  primary  objective  of  this  section  is  to  summarize  available 
information  about  the  following  aquifer  properties: 

1.  The  geologic  structure  of  the  aquifer  and  its  surroundings 

2.  The  extent  of  subsurface  leakage  into  or  out  of  the  aquifer 

3.  The  location  and  significance  of  areas  where  groundwater  (lows  from  the  surface 
into  the  aquifer  (recharge)  or  from  the  aquifer  to  the  surface  (discharge). 

4.  The  distribution  of  hydraulic  head  throughout  the  aquifer 

5.  The  range  of  likely  values  for  aquifer  transmissivity  and  storage  coefficient. 

This  is  the  basic  information  needed  to  construct  flow  nets,  water  budgets,  and 
other  more  complex  models  of  groundwater  flow  in  the  aquifer. 

The  geologic  strata  in  the  San  Juan  basin  can  be  divided,  for  present  purposes, 
into  the  following  groups  (see  Figure  2-3  for  a  typical  geological  cross  section): 

1.  Older  strata  bounded  on  the  top  by  the  Yeso  formation  and  extending  down 
to  Precambrian  rock. 

2.  Strata  of  intermediate  age,  including  the  San  Andres  limestone,  the  Glorieta 
sandstone,  and  the  Chinle  formation. 

3.  Newer  strata  composed  of  alternating  layers  of  shale  and  sandstone  extending 
up  to  recent  alluvium. 

The  uplift  of  the  Zuni  mountains  has  resulted  in’  the  progressive  exposure  of  each 
of  these  strata  along  an  axis  extending  southwest  to  northeast  perpendicular  to  the 
range. 

The  San  Andres  and  Glorieta  formations  are  the  major  sources  of  groundwater 
in  the  region — they  form  a  single  hydrologically  connected  unit  generally  referred 
to  as  the  San  Andres- Glorieta  aquifer.  The  Glorieta  sandstone  has  the  largest  areal 
extent  of  the  two  formations  (see  Figure  2-1).  It  is  exposed  along  the  northern  flank 
of  the  Zuni  mountains  and  gradually  pinches  out  about  80  miles  north.  The  San 
Andres  consists  of  two  limestone  beds  separated  by  an  intervening  sandstone  layer 
similar  in  composition  to  the  Glorieta.  The  San  Andres  outcrops  along  portions 
of  the  northern  flank  of  the  Zunis — in  these  areas  it  is  eroded  and  has  a  relatively 
high  transmissivity.  The  San  Andres  has  a  maximum  thickness  of  about  200  feet, 
as  compared  to  300  feet  for  the  Glorieta,  and  disappears  altogether  in  some  areas 
north  of  the  Zuni  mountains.  It  is  worth  noting  that,  although  the  northern  and 
southern  limits  of  the  San  Andres- Glorieta  aquifer  are  relatively  well-defined,  it’s 
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FIGURE  2-3 

TYPICAL  GEOLOGICAL  CROSS  SECTION  THROUGH  THE  SAN  JUAN  BASIN 


eastern  and  western  boundaries  are  less  clear.  According  to  information  presented 
by  Baars  and  Stephenson  (1977),  these  boundaries  appear  to  lie  well  beyond  the 
area  covered  by  Figure  2-1. 

Available  evidence  on  the  San  Juan  basin  suggests  that  the  San  Andres-  Glorieta 
is  confined,  except  near  outcrop  areas,  by  the  Chinle  shale  above  and  the  siltstone 
and  claystone  beds  of  the  Yeso  formation  below.  The  Chinle  is  a  large  formation 
found  throughout  most  of  northern  Arizona  and  New  Mexico.  The  areal  extent  of 
the  Yeso  appears  to  be  about  the  same  as  the  Glorieta.  Although  the  San  Andres- 
Glorieta  is  surrounded  by  less  permeable  geologic  formations,  there  is  probably 
some  subsurface  movement  of  water  (leakage)  into  or  out  of  the  aquifer.  Both  the 
Chinle  and  Yeso  formations  include  sandstone  beds  which  are  sufficiently  permeable 
to  supply  at  least  modest  amounts  of  water. 

The  direction  and  magnitude  of  leakage  between  the  San  Andres- Glorieta  and 


adjoining  formations  depends  on  the  direction  and  magnitude  of  the  head  gradient 
across  the  formation  boundaries  as  well  as  on  the  Icakancc  value  at  the  boundary. 
Since  both  the  head  gradient  and  leakage  vary  over  the  aquifer  it  is  difficult  to 
make  general  statements  about  leakage.  Investigators  appear  to  disagree  about  the 
direction  of  vertical  gradients  between  the  San  Andrcs-Glorieta  and  its  neighboring 
strata.  Assumed  or  estimated  ieakance  values  vary  from  0.0  (IIEC,  1982)  to  1.0  X 
10  11  sec-1.  (Gcotrans,  1982).  An  approximate  estimate  of  net  leakage  flux  can  be 
obtained  by  computing  a  water  budget  for  the  aquifer.  The  budget  presented  at 
the  end  of  this  section  suggests  that  the  net  direction  of  leakage  is  out  of  the  San 
Andres-Gloricta  into  adjoining  formations. 

Surface  recharge  to  the  San  Andrcs-Glorieta  comes  mostly  from  areas  where 
the  limestone  or  sandstone  layers  of  the  aquifer  are  exposed  or  where  these  layers  are 
overlain  by  permeable  beds  of  the  Chinle  formation,  The  major  recharge  area  lies 
along  the  higher  elevations  of  the  Zuni  mountains,  in  the  shaded  region  indicated 
on  Figure  2-1.  Note  that  the  aquifer  is  unconlined  (has  a  free  surface  for  an 
upper  boundary)  where  recharge  takes  place.  This  implies  that  the  hydraulic  head 
response  near  the  Zuni  mountain  recharge  area  will  be  slower  and  more  damped 
than  in  interior  regions  of  the  aquifer  where  artesian  conditions  prevail. 

Recharge  in  the  Zuni  mountain  outcrop  region  comes  mostly  from  snowmelt 
and  does  not  appear  to  vary  greatly  from  year-to-year.  Most  investigators  seem 
to  accept  Shomaker’s  (1971)  estimate  that  the  long-term  average  (approximately 
steady-state)  annual  infiltration  rate  in  the  recharge  zone  is  about  1.0  inch/year. 
It  should  be  pointed  out  that  this  value  is  speculative  since  supporting  field  ob¬ 
servations  are  limited  (Shomaker,  1971;  Geotrans,  1982).  If  the  1.0  inches/year 
infiltration  rate  is  applied  uniformly  over  a  recharge  area  of  150  square  miles  (a 
conservatively  low  value)  the  estimate  obtained  for  total  annual  recharge  is  8040 
acre-feet/year.  Head-dependent  leakage  from  Bluewater  Lake  and  other  surface 
waters  could  contribute  additional  recharge,  although  these  sources  are  probably 
small. 

Groundwater  discharge  from  the  San  Andres- Glorieta  is  even  more  difficult 
to  estimate  than  recharge,  partly  because  it  is  less  localized  and  varies  more  over 
time.  Temporal  variations  in  total  annual  discharge  from  the  San  Andres- Glorieta 
are  illustrated  by  Shomaker’s  (1971)  estimates  for  the  period  1929-1969  (see  Figure 
2-4).  These  estimates  clearly  demonstrate  the  increase  in  groundwater  usage  which 
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has  occurred  over  the  last  few  decades. 

The  most  convenient  way  to  summarize  available  discharge  information  is 
to  distinguish  two  discharge  categories  -  (1)  small  wells  and  springs  with  nearly 
constant  discharges  and  (2)  large  wells  which  account  for  the  observed  increase 
in  regional  pumpage.  The  first  category  includes  about  90  stock  wells  scattered 
throughout  the  aquifer  and  a  few  springs  on  the  northern  flank  of  the  Zuni  moun¬ 
tains.  Each  of  these  small  sources  discharges  an  average  of  about  1.5  acre-feet/year 
(Shomaker,  1971).  Table  2-1  summarizes  estimated  discharges  for  the  second  cate¬ 
gory  for  two  years,  1968  and  1982.  The  values  given  are  rough  estimates  intended 
for  comparative  purposes  only  and  do  not  necessarily  represent  the  water  rights 
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claimed  by  the  water  users  listed  (claimed  rights  are  generally  larger  than  historical 
pumpage  values).  Table  2-1  shows  that  the  larger  wells  in  Category  2  are  responsible 
for  most  of  the  discharge  from  the  aquifer.  These  wells  tend  to  be  concentrated  in 
a  narrow  strip  running  along  Interstate  40  between  Grants  and  Gallup. 

There  are  relatively  few  surveys  of  the  regional  hydraulic  head  distribution 
in  the  San  Andres-Glorieta  aquifer.  Shomaker  (1971)  conducted  a  study  of  well 
observations  in  the  Fort  Wingate  vicinity  (east  of  Gallup)  which  shows  declines  of 
from  40  to  150  feet  over  the  period  1958  to  1968.  Geohydrology  (1982)  extrapolated 
Shomaker’s  1968  data  several  miles  east  using  more  recent  information  obtained 
from  scattered  well  observations.  The  extrapolated  distribution  is  shown  in  Figure 
2-5  (in  blue)  together  with  some  1958  contours  compiled  by  Gordon  (1961)  for  the 
Grants- Bluewater  area  (in  red).  Gcotrans  (1982)  and  HEC  (1982)  assumed  that  the 
1968  heads  represent  an  approximate  steady-state.  This  assumption  is  supported 
by  data  which  indicate  that  annual  aquifer  discharge  remained  nearly  constant 
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TABLE  2-1 


ESTIMATED  HISTORICAL  DISCHARGES  FROM  THE 
SAN  ANDRES-GLORIETA  AQUIFER 


(all  values  in  acrc-fect  / year) 


CATEGORY  I  Small  Wells  and  Springs 

90  sources  @1.5  acre-fcet 

135 

Subtotal 

135 

CATEGORY  II-Large  Wells  and  Well  Fields 

RANGE 

Annual  Pumpage 
1968  1982 

El  Paso  Natural  Gas 

13 

13 

70 

80 

Independent  users 

13 

13 

120 

220 

Bureau  of  Indian  Affairs 

13 

14 

50 

90 

Transwestern  Pipeline 

14 

13 

75 

75 

Town  of  Thoreau 

14 

13 

40 

75 

Independent  users 

14 

14 

90 

200 

Whispering  Cedars 

14 

15 

50 

200 

Independent  users 

14 

15 

135 

250 

Shell  Oil  Refinery 

15 

15 

600 

600 

El  Paso  Natural  Gas 

15 

16 

200 

225 

El  Paso  Natural  Gas 

15 

17 

200 

275 

U.S.  Army  (Ft.  Wingate) 

15 

17 

— 

200 

Subtotal 

1630 

2490 

AQUIFER  TOTALS 

1765 

2625 

‘ts 


for  several  years  prior  to  1968  (see  Figure  2-4).  Although  the  1968  head  data  do 
not  necessarily  reflect  current  water  levels  in  the  San  Andrcs-Glorieta,  they  do  give 
a  good  qualitative  picture  of  groundwater  movement.  Figure  2-5  clearly  shows  that 
the  general  direction  of  flow  is  northward,  away  from  the  recharge  area,  with  a 
gradient  of  approximately  60  feet/mile  at  the  northern  edge  of  the  aquifer  outcrop. 

Practical  estimates  of  aquifer  transmissivity  and  storage  coefficient  are  usually 
based  on  several  different  sources  of  data.  The  following  paragraphs  briefly  review 
some  relevant  geologic  and  hydrologic  information  available  for  the  San  Andres- 
Glorieta.  An  analysis  of  the  way  this  information  may  be  used  to  develop  model 
inputs  for  the  case  study  is  provided  in  Chapter  4. 

Transmissivity  and  storage  coefficient  are  vertically  averaged  aquifer  parameters 
which  may  be  obtained  by  multiplying  the  hydraulic  conductivity  and  specific 
storage,  respectively,  by  the  groundwater  flow  depth  (see  Section  3.2.2).  Under  ar¬ 
tesian  conditions,  this  depth  is  equal  to  the  aquifer  thickness.  Since  the  San  Andres- 
Glorieta  is  composed  of  two  geological  formations  which  have  different  hydrologic 
properties  it  is  best  to  consider  the  thickness  of  each  formation  when  estimating 
vertically  averaged  parameters.  A.quifer  thickness  maps  computed  by  Gcotrans 
(1982)  indicate  that  the  Glorieta  sandstone  gradually  decreases  in  thickness  from 
the  Zuni  mountains  north  with  a  moderate  bulge  occurring  east  of  Crown  Point. 
The  thickness  of  the  San  Andres  limestone  varies  significantly,  primarily  as  a  result 
of  the  effects  of  erosion  and  fracturing.  Available  well  logs  suggest  that  the  San 
Andres  may  vanish  altogether  in  some  areas  near  the  Plains  well  field.  The  major 
contribution  of  this  formation  probably  occurs  near  the  Zuni  mountain  recharge 
area. 

There  appear  to  b'-  significant  regional  variations  in  the  hydraulic  conductivity 
and  specific  storage  of  the  San  Andres  and  Glorieta  formations.  The  porosity  and 
conductivity  of  the  San  Andres  increases  where  there  has  been  erosion,  fracturing 
and  minor  faulting,  as  has  been  observed  in  the  Grants-Bluewater  region.  The 
Glorieta  sandstone  is  generally  less  heterogeneous  than  the  San  Andres  formation, 
although  its  conductivity  and  porosity  apparently  increase  where  the  limestone  has 
eroded.  This  tends  to  compensate  somewhat  for  the  disappearance  of  the  more 
conductive  San  Andres.  Some  faults  and  anomalies  found  near  the  Zuni  mountains 
may  block  groundwater  flow,  particularly  in  the  direction  parallel  to  the  axis  of  the 
Zuni  range.  Although  these  structural  features  could  influence  the  results  of  local 
pump  tests  they  probably  have  little  effect  on  regional  (aquifer  scale)  flow  patterns, 
at  least  north  of  the  Zuni  outcrop. 

When  the  above  factors  are  considered  together  a  qualitative  description  of 
aquifer  parameter  variations  emerges.  Transmissivity  is  highest  in  the  southern  por¬ 
tion  of  the  aquifer,  particularly  near  Grants,  and  gradually  decreases  to  the  north. 
A  similar  but  less  dramatic  variation  probably  applies  to  the  storage  coefficient. 
Locally  high  values  of  transmissivity  may  be  observed  where  the  San  Andres  is 
present  if  fracturing  or  solution  porosity  are  significant.  Finally,  aquifer  properties 
are  relatively  homogeneous  in  areas  where  the  San  Andres  is  absent. 
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Figure  2-5 

MEASURED  HEAD  CONTOURS  IN  THE  SAN  ANDRES- GLORIETA  AQUIFER 


Although  this  qualitative  description  is  useful,  it  does  not  provide  actual  values 
for  the  aquifer  parameters.  These  values  must  be  estimated  using  one  of  the 
procedures  mentioned  in  Chapter  3.  For  an  aquifer  the  size  and  complexity  of 
the  San  Andres-Glorieta  the  most  practical  estimation  techniques  are  pump  test 
analysis  and  regional  model  calibration.  These  methods  are  not  mutually  exclusive 
but  may  be  combined  to  make  best  use  of  available  information. 

Pump  test  results  for  the  San  Andres-  Glorieta  are  available  from  Shomaker’s 
(1971)  report  on  the  Ft.  Wingate  vicinity  and  from  Geohydrology  (1982).  The 
Shomaker  results  give  transmissivities  ranging  from  5  to  3700  ft2/day,  depending 
on  location.  The  storage  coefficients  (based  only  on  a  few  wells)  range  between 
3.0  X  10~5  and  1.3  X  10~4  (unitless).  It  should  be  noted  that  these  estimates  come 
from  a  variety  of  sources  and  are  often  based  on  short-term  (i.e.,  limited  range)  tests 
on  wells  which  may  be  partially  or  multiply  completed.  They  are  useful  primarily 
for  defining  an  approximate  range  of  likely  values. 

The  pump  tests  reported  by  Geohydrology  (1982)  were  conducted  at  the  pro- 


posed  Plains  Electric  well  Held  and  have  therefore  received  considerable  scrutiny  by 
modelers  concerned  with  the  Plains  pumping  proposal.  Transmissivity  estimates 
based  on  these  tests  vary  from  80  ft'/day  (Dames  and  Moore,  1982)  to  3000  ft“/day 
((.’eohydroiogy,  1982)  although  most  investigators  appear  to  now  agree  on  the 
narrower  range  of  150  to  450  ft“/day.  Storage  coefficient  estimates  vary  from 
1.0  *  10  1  (Dames  and  Moore,  1982)  to  5.0  X  lO-'1  (Geohydrology,  1982). 

The  primary  data  available  for  regional  model  calibration  are  the  1968  head 
contours  shown  in  Figure  2-5.  These  heads  were  used  by  both  Gcotrans  (1982)  and 
MIX’  (1982)  to  guide  adjustments  of  transmissivity  and  other  model  inputs  such  as 
recharge  and  leakage.  Geohydrology  (1982)  used  head  observations  from  the  area 
near  the  Shell  refinery  well  Geld  for  a  similar  purpose.  While  the  Shell  observations 
are  not  regional  in  nature,  they  do  provide  useful  information  about  the  long-term 
response  to  stress  in  an  important  part  of  the  study  area. 

Geohydrologic  information  presented  in  this  section  is  conveniently  summarized 
with  a  simplified  water  budget  for  the  San  Andres-Glorieta  aquifer.  The  major  com¬ 
ponents  of  the  water  budget  are  recharge  (supply),  discharge  and  leakage  (demand), 
and  storage  change.  The  recharge,  discharge,  and  head  data  described  earlier 
provide  the  information  needed  to  estimate  two  water  budgets  -  a  quasi-steady- 
state  budget  for  1968  and  a  non-steady-state  budget  for  1982.  Each  of  these  is 
summarized  in  Table  2-2. 

The  recharge  is  based  on  the  1.0  inch/year  infiltration  rate  suggested  by 
Shomakcr  (1971)  and  is  assumed  to  be  relatively  constant.  Recharge  from  Blue- 
water  Lake  and  other  surface  waters  (if  any)  is  not  included.  Discharge  values 
are  obtained  from  Table  2-1.  Leakage  is  selected  as  the  unknown  element  of  the 
1968  steady-state  water  budget  since  it  is  the  most  difficult  to  estimate.  If  the 
1968  leakage  is  adjusted  to  give  a  water  balance  the  result  is  6475  acre-feet  /year. 
This  is  an  estimate  of  the  net  amount  of  subsurface  water  leaving  the  San  Andres- 
Glorieta  aquifer  during  1968.  A  water  balance  can  be  obtained  for  1982  if  the  1968 
leakage  value  is  assumed  to  apply  and  storage  change  is  selected  as  the  unknown. 
This  gives  an  estimated  1982  storage  decrease  of  860  acre-feet.  This  decrease  is  due 
both  to  compressibility  (artesian)  effects  and  to  water  table  declines  in  unconfined 
portions  of  the  aquifer. 

Although  the  estimates  given  in  Table  2-2  are  speculative,  they  clearly  suggest 
that  leakage  is  an  important  component  of  the  water  budget.  The  estimated 
recharge  would  have  to  be  decreased  significantly  (by  a  factor  of  five)  to  change 
this  conclusion.  Such  a  decrease  could  be  obtained  only  by  limiting  infiltration 
to  unreasonably  low  values  or  by  greatly  reducing  the  generally  accepted  area  of 
the  recharge  region.  A  rough  check  on  the  recharge  estimate  can  be  obtained 
by  using  the  head  gradient  of  60  feet/mile  observed  at  the  northern  edge  of  the 
Zuni  mountain  outcrop.  If  this  value  is  multiplied  by  the  length  of  the  recharge 
area  (approximately  50  miles)  and  a  conservatively  low  transmissivity  value  of  300 
ft2/day,  the  resulting  estimate  of  subsurface  flux  entering  the  aquifer  from  the  Zuni 
outcrop  is  7500  acre-feet/year.  This  is  consistent  with  Table  2-2. 
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TABLE  2  2 


SIMPLIFIED  WATER  BUDGETS  FOR  THE 
SAN  A N 1) R ES- G LOR I ETA  AQUIFER 

(all  values  in  acre-fcet  /  year) 

SUPPLY  (into  aquifer) 

1968 

1982 

(  water  levels')  f 

water  levels 

l  steady  J  1 

falling 

Recharge  from  Zuni  mountain  outcrop 

8040 

8040 

Subtotals 

8040 

8040 

DEMAND  (out  of  aquifer) 

Discharge  from  wells  and  springs 

1765 

2625 

Leakage  to  adjoining  aquifers 

6275 

6275 

Subtotals 

8040 

8900 

STORAGE  CHANGE 
Annual  decrease  in  storage 


860 


3.  BASIC  MODELING  CONCEPTS 


The  case  study  described  in  the  preceding  chapter  has  a  number  of  features 
which  are  commonly  encountered  in  ground  water  management.  The  basic  prob¬ 
lem  is  to  predict  impacts  relatively  far  into  the  future  (forty  years)  using  a  limited 
amount  of  ambiguous  data.  The  economic  stakes  are  large,  both  for  the  power  com¬ 
pany  and  for  local  residents  who  may  be  adversely  affected  by  increased  groundwater 
pumpage.  Mathematical  modeling  seems  to  be  a  reasonable  approach  to  use,  al¬ 
though  the  limitations  of  the  available  data  base  cast  some  doubt  on  the  accuracy 
of  the  any  long-term  prediction. 

This  chapter  considers  in  detail  how  groundwater  models  may  be  applied  in 
case  studies  such  as  the  San  Andres-Glorieta.  Individual  sections  describe  the 
various  tasks  which  are  required  to  derive  credible  predictions  from  uncertain 
field  data.  These  tasks  were  carried  out  in  all  of  the  modeling  studies  examined 
in  Chapter  4,  although  the  specific  methods  used  and  answers  obtained  differed 
considerably.  As  the  discussion  moves  from  general  concepts  to  specific  criticisms 
it  should  become  apparent  how  the  methods  outlined  here  could  be  used  to  evaluate 
other  modeling  studies  directed  toward  other  management  problems. 

A  typical  modeling  study  includes  three  basic  elements  -  definition  of  the 
problem  to  be  investigated,  formulation  of  a  modeling  strategy  for  solving  this 
problem,  and  actual  application  of  the  model.  These  can  be  divided  into  more 
specific  tasks  as  follows: 


Problem  Definition 

1.  Identify  the  management  problem  to  be  investigated  (e.g.,  evaluation  of  the  im¬ 
pacts  of  aquifer  development)  and  formulate  a  set  of  preliminary  specifications 
for  the  modeling  study. 

2.  Review  existing  information  on  the  site  (e.g.,  geological  reports,  well  data,  etc.). 
Develop  a  feeling  for  the  quantity  and  quality  of  the  data  available.  Construct 
a  qualitative  description  of  groundwater  flow  in  the  aquifer  of  interest. 


Model  Formulation 

1.  Develop  a  detailed  mathematical  description  of  the  aquifer  system.  This  descrip¬ 
tion  should  account  for  the  physical  factors  most  relevant  to  the  problem  of 
interest  and  should  be  specific  enough  to  indicate  the  type  of  computer  program 
required. 


Model  Application 

1.  Select  a  computer  program  for  solving  the  modeling  problem. 

2.  Determine  the  level  of  spatial  and  temporal  discretization  required  in  the 
model’s  inputs  and  outputs  and  set  up  a  simulation  network. 

3.  Estimate  all  model  inputs  from  available  data  sources. 

4.  Check  the  model’s  prediction  accuracy  using  comparisons  of  model  predictions 
with  field  data  and/or  sensitivity  analyses. 

5.  Use  the  model  to  investigate  the  management  problem  posed  at  the  beginning 
of  the  study. 

The  discussion  of  modeling  concepts  presented  in  this  chapter  follows  the  organiza¬ 
tion  outlined  above.  Section  3.1  defines  the  impact  evaluation  problem  and  suggests 
how  a  qualitative  description  of  groundwater  flow  can  be  constructed  from  com¬ 
monly  available  data  sources.  Section  3.2  reviews  the  general  principles  used  to 
develop  a  mathematical  groundwater  model.  The  modeler  should  have  a  working 
knowledge  of  these  principles  even  if  he  plans  to  use  an  “off-the-shelf”  computer 
program  since  many  practical  modeling  decisions  depend  on  the  properties  of  the 
governing  equations.  Section  3.3  examines  a  number  of  application-oriented  issues, 
including  the  important  problems  of  input  estimation  and  accuracy  evaluation. 
These  three  sections  provide  most  of  the  technical  background  needed  for  the  case 
study  analysis  of  Chapter  4. 


3.1  PROBLEM  DEFINITION 

The  management  issue  addressed  in  the  San  Andres- Glorieta  case  study  is  the 
problem  of  evaluating  the  impact  of  pumpage  from  a  limited  groundwater  resource. 
This  problem  arises  frequently  in  water  supply  planning,  particularly  in  arid  regions 
such  as  the  San  Juan  basin,  where  groundwater  is  the  major  supply  source.  A 
comprehensive  impact  evaluation  should  probably  include  some  consideration  of 
the  economic,  social,  water  quality,  and  hydrologic  consequences  of  development. 
In  practice,  attention  is  usually  focused  on  those  impacts  which  can  be  most  readily 
measured  -  changes  in  aquifer  water  levels  and  changes  in  groundwater  quality. 
This  report  considers  only  hydrologic  impacts,  primarily  because  of  a  need  to  limit 
the  scope  of  the  presentation.  It  should  be  noted,  however,  that  most  of  the  general 
concepts  discussed  here  are  also  relevant  to  investigations  of  water  quality. 

When  a  groundwater  hydrologist  decides  to  use  a  mathematical  model  for 
impact  analysis,  he  should  identify  as  clearly  as  possible  the  type  of  information 
he  expects  to  obtain.  One  way  to  do  this  is  to  attempt  to  answer  the  following 
questions: 

1.  What  variables  should  be  used  to  measure  hydrologic  impact?  Depth- to- water, 
change  in  storage,  drawdown,  or  other  less  common  indices? 
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2.  Where  lire  the  boundaries  of  the  aquifer?  Do  they  extend  beyond  the  region  of 
most  interest  for  impact  evaluation? 

3.  Where  are  predictions  required?  How  much  spatial  resolution  is  needed?  Do 
some  areas  warrant  a  more  detailed  analysis  than  others? 

4.  How  far  in  the  future  arc  predictions  required?  How  much  temporal  detail  is 
needed?  Are  certain  seasons  or  years  particularly  important? 

5.  How  accurate  must  the  model’s  predictions  be  in  order  to  be  useful  for  decision¬ 
making? 

The  answers  to  these  questions  constitute  a  set  of  preliminary  specifications  for  the 
modeling  study. 

It  is  generally  helpful  to  review  available  aquifer  data  before  a  mathematical 
model  is  formulated.  In  the  San  Andres-Glorieta  case  study,  useful  sources  of 
information  included  geological  reports  and  maps  prepared  by  the  U.S.  Geological 
Survey,  consultant  or  government  reports  documenting  earlier  studies,  and  well  logs 
or  water  level  observations  compiled  locally.  A  preliminary  “desk-top”  analysis  of 
some  of  this  information  provides  a  good  feeling  for  regional  groundwater  behavior 
and  may  even  indicate  that  a  detailed  modeling  investigation  is  not  necessary. 

A  good  example  of  a  valuable  desk- top  calculation  is  the  San  Andres-Glorieta 
aquifer  water  budget  presented  in  Table  2-2.  This  type  of  simple  water  budget 
quickly  identifies  serious  data  gaps  which  must  be  filled  before  a  complete  modeling 
effort  can  be  undertaken.  If  water  level  data  are  available  at  two  or  more  times 
they  can  be  used  to  check  the  storage  change  computed  by  subtracting  supply  and 
demand  figures.  This  can  reveal  inconsistencies  in  data  and  assumptions  which  are 
best  discovered  early  in  the  study. 

There  are  a  number  of  other  simple  analyses  that  can  be  performed  prior  to 
modeling.  These  include  construction  of  flow  nets  and  calculations  of  well  field 
drawdowns  using  Theis  curves  or  other  graphical  methods  (see  Freeze  and  Cherry, 
1979,  and  Chapter  4  of  this  report  for  some  examples).  The  results  of  several 
different  desk-top  computations  can  be  combined  to  give  an  overall  picture  of 
groundwater  flow  in  the  aquifer  of  interest.  This  may  later  be  used  to  check  the 
credibility  of  the  model’s  predictions. 


17 


3.2  MODEL  FORMULATION 

3.2.1  Basic  Concepts  of  Groundwater  Flow 

An  aquifer-scale  water  budget  provides  a  useful  but  highly  aggregated  descrip¬ 
tion  of  subsurface  Dow.  A  more  detailed  analysis  is  needed  to  predict  local  changes 
in  water  levels  and  to  account  for  spatial  heterogeneity  in  applications  such  as  the 
San  Andrcs-Clorieta  case  study.  A  detailed  analysis  requires  the  water  budget  ap¬ 
proach  to  be  applied  individually  to  many  small  pieces  (or  elements)  of  the  aquifer. 
The  smallest  aquifer  element  which  can  be  analyzed  with  traditional  groundwater 
[low  theory  is  the  so-called  “representative  elementary  volume”,  commonly  ab¬ 
breviated  REV  (see  Freeze  and  Cherry,  1979).  This  can  be  thought  of  as  a  volume 
element  whose  dimensions  are  small  compared  to  the  entire  aquifer  but  large  com¬ 
pared  to  individual  pores  in  the  aquifer  soil  matrix.  The  actual  size  of  an  elementary 
volume  depends  on  the  type  of  aquifer — typical  dimensions  for  the  San  Andres- 
Glorieta  vary  from  a  few  inches  for  a  fine  sandstone  to  hundreds  of  feet  for  a 
fractured  limestone.  It  is  mathematically  convenient  to  derive  water  balances  and 
other  equations  describing  groundwater  flow  at  the  REV  scale  and  then  to  integrate 
these  equations  to  obtain  solutions  for  larger  regions. 

Some  of  the  basic  definitions  used  in  groundwater  flow  derivations  are  il¬ 
lustrated  in  Figure  3-1.  This  figure  shows  a  hypothetical  experimental  apparatus 
designed  to  measure  the  properties  of  a  saturated  porous  medium  (e.g.,  sand).  The 
sand  is  confined  in  a  cylinder  whose  dimensions  are  large  compared  to  the  repre¬ 
sentative  elementary  volume  illustrated  in  the  detail.  Water  flows  from  an  inlet 
tank  with  a  surface  elevation  hj  through  the  sand  into  an  outlet  tank  with  a  lower 
surface  elevation  /12.  The  respective  inflow  and  outflow  rates  (in  units  of  l3/t)  are 
designated  Qin  and  Qout-1 

The  water  balance  for  the  cylinder  simply  states  that  the  difference  between 
the  inflow  and  outflow  at  any  time  is  equal  to  the  rate  of  change  in  the  mass  of 
water  stored  in  the  sand.  This  may  be  expressed  mathematically  in  the  following 
way: 

(3-D 

where  p  is  the  density  of  water  (with  units  m/l 3),  n  is  the  unitless  porosity  defined 
in  Figure  3-1,  and  V  is  the  volume  of  the  cylinder.  The  product  pnV  is  the  total 
mass  of  water  in  the  cylinder  at  a  given  time. 

'In  this  discussion,  the  following  symbols  arc  used  to  define  units: 
l  —  length 
m  =  mass 
t  =  time 
/  =  forcc(  ml/t1) 
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FIGURE  3-1 

EXPERIMENTAL  APPARATUS  USED  TO  ILLUSTRATE  BASIC  CONCEPTS 

OF  GROUNDWATER  FLOW 


At  the  smaller  REV  scale  an  equivalent  water  balance  equation  may  be  written 
in  differential  form: 

d{pn)  =  d[pq{l,t)\ 
dt  di 

Here  q(l,  t )  is  defined  as  the  specific  discharge,  or  flow  per  unit  area,  through  the 
REV  at  location  l  and  time  t  (see  Figure  3-1).  Partial  derivatives  are  used  because 
the  discharge  depends  on  more  than  one  independent  variable.  Note  that  the  minus 
sign  is  required  to  insure  that  the  mass  decreases  if  the  flow  increases  in  the  direction 
of  positive  l. 

The  dependence  of  p  and  n  on  pressure  are  accounted  for  reasonably  well  by 
the  following  expression  (Freese  and  Cherry,  1979): 


(3-2) 


I 

I. 


where  h  is  the  hydraulic  head  (with  units  /)  and  S,,  is  a  soil  parameter  called  the 
spccilic  storage  (with  units  / ' Equations  (3-2)  and  (3-3)  may  be  combined  to  give 
a  compact  expression  for  the  dynamic  mass  balance  at  the  REV  scale: 

=  (3-4) 

at  dt  K  ’ 

This  is  sometimes  called  the  mass  continuity  equation. 

The  derivation  leading  to  Equation  (3-3)  provides  a  convenient  expression  for 
the  specific  storage: 

S„  =  pg(a  +  n/3)  (3-5) 

Here  pg  is  the  specific  gravity  of  water  (with  units  ///3),  g  is  the  acceleration  of 
gravity  (with  units  l/t2),  and  a  and  /?  arc  the  aquifer  and  fluid  compressibility  coeffi¬ 
cients  (with  units  l2/ f).  Spatial  variations  in  specific  storage  are  due  primarily  to 
variations  in  aquifer  porosity  and  compressibility;  the  other  parameters  in  Equation 
3-5  are  nearly  constant. 

The  unknowns  in  the  continuity  equation  are  related  by  Darcy’s  law,  an  em¬ 
pirical  principle  which  states  that  specific  discharge  is  directly  proportional  to  the 
head  gradient: 


ah(i,t) 


(3-6) 


Here  K  represents  the  hydraulic  conductivity  of  the  soil  (with  units  f/t).  This 
parameter  depends  on  the  structure  of  the  soil  or  rock  matrix  and  may  vary  over 
many  orders  of  magnitude.  It  can  be  estimated  from  the  apparatus  shown  in  Figure 
3-1  if  the  specific  discharge  and  hydraulic  head  are  held  constant. 

Darcy’s  law  allows  the  REV  mass  balance  to  be  written  as  a  single  differential 
equation  in  one  unknown,  the  hydraulic  head: 


£MU)_£  afc[u) 

s'  at  arK  ai  1 


(3-7) 


This  expression  may  be  generalized  to  three-dimensional  situations  to  give  the 
fundamental  equation  of  groundwater  flow  in  various  coordinate  systems.  The 
well-known  cartesian  version  of  this  equation  obtained  for  isotropic  (direction- 
independent)  hydraulic  conductivity  and  specific  storage  may  be  written: 


C  _  w 

~  ^  r. 


*[  Kd± 

dy  dy 


.±[ Kd± 

dz  dz 


(3-8) 


Here  it  should  be  understood  that  the  head  depends  on  the  coordinates  of  the  REV 
(z,  y,  and  z )  as  well  as  on  time.  In  saturated  soils  the  isotropic  parameters  St  and 
K  are  functions  of  location  only. 


•  .*•  **  ^  .**  . 


3.2.2  The  Vertically  Averaged  Flow  Equation 


Many  aquifers  used  for  water  supply  arc  thin  geological  strata  which  extend 
horizontally  for  many  miles  but  may  be  only  a  few  hundred  feet  thick.  The  San 
Andres-Gloricta  is  a  typical  example,  as  is  illustrated  by  the  geological  cross  section 
shown  in  Figure  2-3.  In  such  aquifers,  variations  in  head  over  the  smaller  vertical 
dimension  may  be  ncgligable  compared  to  variations  over  the  larger  horizontal 
dimension.  In  order  to  consider  some  of  the  practical  implications  of  this,  it 
is  useful  to  distinguish  two  types  of  head  observations.  The  three-  dimensional 
head  h(x,y,z,t)  appearing  in  Equation  (3-8)  can  be  thought  of  as  the  water  level 
measured  in  a  well  located  at  horizontal  coordinates  x  and  y  which  has  a  single 
very  short  well  screen  installed  at  vertical  elevation  z.  If  the  screen  elevation  were 
changed  to  z' ,  the  water  level  in  the  well  could  also  change,  although  the  difference 
may  be  ncgligable.  If  the  well  screen  extended  throughout  the  entire  depth  of  the 
aquifer,  from  a  lower  elevation  of  Z^(x, y)  to  an  upper  elevation  of  Zu(x,y),  the 
water  level  observed  in  the  well  would  be  an  average  of  the  heads  at  all  elevations 
from  Zi  to  Z\j.  This  vertically  averaged  head  (written  h(x,y,t ))  may  be  defined 
mathematically  as: 

h(x ,  y,  t)  =  — - -r  /  h(x,  y,  z,  t )  dz 

D(x,  y,  t)  J Zl(z,v) 

where  D[x,y,t)  =  Zy[x,y,t)  —  Zj£xKy)  is  the  saturated  flow  depth.  Note  that  Zy 
and  D  will  depend  on  time  if  the  aquifer  is  unconfined  (i.e.  if  the  upper  boundary 
surface  is  free  to  rise  or  fall).  If  the  aquifer  is  confined,  Z\j  and  D  will  be  constant 
at  any  given  location. 

If  vertical  variations  throughout  the  aquifer  depth  are  small  compared  to 
horizontal  variations  across  the  aquifer,  then  the  unaveraged  head  h(x,y,z,t )  will 
be  close  to  the  averaged  head  Ti(x,  y,  t }  at  any  elevation  z.  In  this  case,  there  is 
really  no  need  to  solve  the  three-dimensional  flow  equation  to  obtain  the  unaveraged 
head.  Instead,  we  can  solve  a  simpler  two-dimensional  equation  which  gives  the 
vertically  averaged  head  at  any  horizontal  location  (x,  y).  This  vertically  averaged 
head  is  effectively  treated  as  an  approximation  to  the  unaveraged  head  at  any  depth 
between  Zi  and  Z\j. 

The  relationship  between  averaged  and  unaveraged  head  is  illustrated  in  Figure 
3-2  for  a  simple  wedge-shaped  aquifer.  The  upper  three-dimensional  portion  of  the 
figure  shows  the  five  and  ten  foot  contours  plotted  as  functions  of  all  three  spatial 
coordinates.  Note  that  the  head  decreases  slightly  as  z  increases  (for  constant  x  and 
y),  suggesting  a  small  upward  flow  (from  high  to  low  head).  The  vertical  variation 
in  head  is  clearly  small  compared  to  the  large  horizontal  variation  observed  in 
the  x  direction.  The  lower  two-dimensional  portion  of  the  figure  shows  the  five 
and  ten  foot  contours  obtained  by  integrating  the  three-dimensional  head  over 
the  z  coordinate.  This  two-dimensional  description  captures  the  major  horizontal 
variations  which  are  of  primary  concern  in  a  regional  impact  evaluation. 
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The  two-dimensional  How  equation  which  gives  the  vertically  averaged  head 
may  be  obtained  by  integrating  Equation  (3-8)  over  the  i  coordinate  (Pindcr  and 
Cray,  1977).  The  result  is  the  following  equation,  which  applies  at  any  location 
(x,y)  in  the  horizontal  plane: 


^ dh  d  d  „,dh 

S  -tt:  =  w~T  —  +  jr~T  —  -f  qu  -  Q,„ 

at  ox  ox  ay  oy 

The  new  variables  appearing  in  this  expression  are  defined  as  follows: 


(3-9) 


h  =vertically  averaged  hydraulic  head(I) 

S  =S.hD  =  storage  coefficient  (unitless) 

T  =KD  =  transmissivity  (l/t) 

D  ^saturated  flow  depth  (/) 

qh  =net  vertical  flux  per  unit  area  crossing  the  upper  and  lower 
boundaries  of  the  three-dimensional  aquifer  (l/t) 

Qw  =wcll  pumpage  per  unit  area  (l/t) 

All  of  the  groundwater  models  discussed  in  this  report  arc  based  on  Equation  (3-9). 

The  vertical  boundary  flux  <p4  accounts  for  all  water  entering  or  leaving  the 
aquifer  across  its  upper  or  lower  boundaries.  This  includes  surface  recharge,  evapo- 
transpiration,  and  leakage,  as  well  as  the  change  in  storage  which  results  if  the 
upper  aquifer  boundary  is  a  free  surface  (i.e.,  if  the  aquifer  is  unconfined).  When 
all  of  these  flux  components  arc  considered  <p,  may  be  written  as  follows  : 

qh  —  qi  +  qu  +  Ql  +  Qf  (3-10) 

where  q,  is  the  net  infiltration  rate  (surface  recharge  minus  evapotranspiration),  qu 
and  qi  are  the  leakage  fluxes  across  the  upper  and  lower  boundaries,  and  qj  is  the 
effective  flux  due  to  changes  in  the  free  surface  at  the  upper  aquifer  boundary. 

Leakage  may  be  treated  as  if  it  were  a  recharge  flux,  or  it  can  be  described  by 
the  following  version  of  Darcy’s  law  (written  for  the  upper  boundary): 


9u  —  Kv 


hu-Ti 


(3-11) 


This  equation  assumes  that  the  leakage  passes  through  a  confining  bed  of  thickness 
Azu  (with  units  /)  and  vertical  hydraulic  conductivity  Kvu  (with  units  l/t).  The 
head  on  the  far  side  of  the  upper  confining  bed  is  hu  (assumed  to  be  known)  and 
storage  in  the  bed  is  neglected.  Sometimes  the  ratio  Kvu/&Zv  (with  units  t~l)  is 
called  the  leakance. 

If  a  portion  of  the  aquifer’s  upper  boundary  is  unconfined,  this  boundary  (called 
the  free  surface  or  water  table)  may  rise  or  fall  in  response  to  pressure  changes. 
When  the  water  table  rises,  groundwater  is  stored  in  the  unsaturated  soil  above. 
Conversely,  groundwater  is  released  from  the  draining  soil  column  when  the  water 
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FIGURE  3-2 

USING  VERTICAL  AVERAGING  TO  OBTAIN  A  TWO-DIMENSIONAL  FLOW  MODEL 


table  falls.  The  amount  of  water  released  (in  mass  per  unit  area)  in  response  to  a 
change  in  the  height  of  the  water  table  is  given  by: 

Am  =  —p[n  —  0r)Az  (3—12) 

where  0r  is  the  specific  retention,  or  fraction  of  water  retained  in  the  unsaturated 
soil  column,  and  A z  is  the  change  in  water  table  elevation.  Since  the  change  in 
elevation  is  approximately  equal  to  the  change  in  the  vertically  averaged  head,  the 
effective  flux  of  water  due  to  water  table  fluctuations  is  given  by: 


„  _  gd* 

V  ~  o7 


(3-13) 


Here  the  term  (n  —  6r)  has  been  replaced  by  Sv,  a  unitless  aquifer  parameter  called 
the  specific  yield.  It  should  be  noted  that  the  San  Andres- Glorieta  aquifer 


is  confined  everywhere  except  in  limited  areas  along  the  southern  outcrop  region. 
Some  simple  methods  for  dealing  with  the  unconfincd  conditions  found  in  this  region 
are  described  in  Chapter  4. 

If  all  the  effects  contributing  to  the  flux  of  water  across  the  aquifer’s  boundaries 
are  included,  the  depth  averaged  flow  equation  may  be  written  as: 

,  .dh  d  -.dh  d  _ dh 

(S  +  s„)^  =  rJ-a-x  +  -  +  UK  -  h)  (3— ,4) 

+  —  h)  +  m  —  Qw 

Here  Ln  and  L\  arc  leakance  coefficients  for  the  upper  and  lower  boundaries,  and 
hu  and  h/  arc  the  heads  in  the  upper  and  lower  confining  layers.”  Note  that  the 
compressibility  and  water  table  storage  change  terms  arc  combined  on  the  left-hand 
side  of  the  equation.  Because  storage  changes  due  to  compressibility  are  minor 
compared  to  those  due  to  water  table  fluctuations,  S  is  usually  several  orders  of 
magnitude  smaller  than  Sv  and  can  be  neglected  where  the  upper  boundary  is 
unconfined.  Sometimes  the  sum  S  +  Sy  is  called  the  storage  coefficient  and  is 
written  as  S.  This  coefficient  is  then  varied  over  a  wide  range  of  values  to  account 
for  transitions  from  confined  to  unconfined  conditions  (this  is  one  of  the  approaches 
taken  in  Chapter  4). 

3.2.3  Integration  of  the  Flow  Equation  over  Space  and  Time 

The  three  dimensional  and  vertically  averaged  groundwater  equations  discussed 
in  the  preceding  sections  apply  at  a  single  location  (the  centroid  of  an  REV)  and  at 
a  single  time.  Since  these  equations  contain  derivatives  of  the  unknown  head  they 
must  be  integrated  over  an  appropriate  spatial  region  and  time  period  before  an 
explicit  head  solution  can  be  obtained.  The  integration  process  effectively  extends 
the  range  of  the  equation  from  a  single  point  to  the  entire  aquifer. 

Whenever  a  differential  equation  is  integrated  one  or  more  arbitrary  constants 
are  introduced.  A  simple  example  is  the  following  equation  giving  the  position  i  of 
a  particle  moving  in  a  straight  line  at  a  fixed  velocity  V : 


dt 


(3-15) 


If  this  equation  is  integrated  over  the  time  period  from  t  =  0  to  t  =  T,  the  result 
is: 

x{t)  =  Vt  +  C  0<t<T  (3-16) 

This  expression  is  a  solution  to  the  original  equation  for  any  constant  value  of 
C — i.e.,  the  particle  can  be  moving  with  velocity  V  at  any  location  x. 

The  non-uniqueness  of  the  solution  is  eliminated  by  specifying  an  auxiliary 
condition  or  constraint  which  defines  the  location  x(0)  of  the  particle  at  the  time 
t  =  0.  Equation  (3-16)  may  then  be  solved  for  C,  giving  the  unique  solution: 

2Thc  overbar  used  earlier  to  denote  vertically  averaged  head  has  been  dropped  here  to  simplify 
notation. 
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x(t)  =  Vt  +  i(0)  0  <t<T  (3-17) 

The  position  x(0)  is  an  initial  condition  for  the  governing  differential  equation. 

Groundwater  simulation  models  perform  numerical  integrations  which  are  similar 
in  general  concept  to  the  above  example.  These  models  can  provide  unique  head 
solutions  only  if  a  sufficient  number  of  auxiliary  conditions  are  specified  by  the 
modeler.  Some  simple  examples  can  help  to  explain  how  these  conditions  arise  in 
practice.  The  example  presented  below  examines  the  use  of  an  initial  condition  in 
a  hypothetical  transient  flow  problem. 


Example:  Initial  Conditions 

Consider  the  rectangular  aquifer  shown  in  Figure  3-3.  This  aquifer  has  a  con¬ 
stant  depth  and  is  isolated  from  its  surroundings  by  flow  barriers  (e.g.,  geological 
faults)  on  all  four  sides  and  by  confining  layers  above  and  below.  A  small  leakage 
flux  moves  uniformly  across  the  upper  boundary  into  or  out  of  an  adjoining  forma¬ 
tion.  The  aquifer  geometry  of  Figure  3-3  suggests  that  the  head  distribution  can 
be  adequately  described  by  the  two-dimensional  vertically  averaged  groundwater 
equation  (Equation  3-14).  The  transmissivity  is  assumed  to  be  a  constant  and  the 
pumpage  and  infiltration  rates  are  assumed  to  be  zero.  The  upper  boundary  leakage 
is  described  by  Equation  (3-11),  with  the  leakance  and  adjoining  aquifer  heads  as¬ 
signed  constant  values.  The  initial  head  is  assumed  to  have  a  uniform  value  ho 
throughout  the  aquifer. 

Since  the  aquifer’s  geometry,  physical  properties,  and  inputs  are  completely 
uniform  in  space  there  is  no  reason  to  expect  the  head  to  vary  over  the  x  and 
y  coordinates,  although  it  may  change  over  time.  It  is  therefore  reasonable  to 
tentatively  assume  that  the  x  and  y  derivatives  in  Equation  (3-14)  are  zero.  This 
assumption  may  later  be  checked  by  substituting  the  resulting  head  distribution 
into  Equation  (3-14)  to  insure  that  it  is  indeed  a  solution. 

If  h  depends  only  on  time,  the  groundwater  equation  becomes  an  ordinary 
differential  equation: 

S~=Lu{hv-h)  (3-18) 

This  may  be  integrated  using  standard  methods  to  give  (see  Hildebrand, 1976): 


h{t)  =  Ce^1  +  hu\l  -  e~^l\  (3-19) 

Here  C  is  an  unknown  constant  of  integration. 

Equation  (3-19)  gives  a  head  distribution  which  satisifes  the  vertically  averaged 
groundwater  equation  for  any  value  of  C  (this  is  easily  checked  by  direct  substitu- 


aV. 


■*  s 


u. 


ymferm  (T  \ 
-tvyisnitssiv/rkj 


t  ,  row  barrier 

a)  plan  view 

unrfcr'm  leakaqe  (q«) 

t  t  x 

A  bfpe^  con-Pln  inq  layen 


\  owe/-  \ayen 


b)  vevHcal  elevation 


FIGURE  3-3 

PLAN  AND  ELEVATION  OF  THE  AQUIFER  ANALYZED  IN  EXAMPLE  1 


tion).  This  non-uniqueness  is  eliminated  by  imposing  the  initial  conditions  denned 
earlier.  If  the  head  at  time  t  =  0  is  equal  to  ho>  then  Equation  (3-19)  at  this  time 
becomes: 

h(0)  =  h0  =  C  (3-20) 

The  integration  constant  is  equal  to  the  initial  head.  Consequently,  the  unique 
solution  to  the  problem  of  Figure  3-3  is: 

h{t)  =  h0e~ +  hu  1  —  «" (3-21) 

This  head  applies  uniformly  throughout  the  aquifer. 

The  solution  of  Equation  (3-18)  implies  that  the  head  will  rise,  fall,  or  remain 
constant,  depending  on  the  relationship  between  ho  and  hu.  These  three  possibilities 
are  illustrated  in  the  plots  of  head  vs.  time  presented  in  Figure  3-4.  Note  that  the 


26 


head  in  every  case  eventually  approaches  the  constant  value  hu.  This  implies  that 
water  will  leak  either  into  or  out  of  the  aquifer  until  the  two  heads  h[t)  and  hu  are 
equal.  As  this  occurs,  the  leakage  gradually  diminishes  and  the  aquifer  approaches 
steady-state  (i.e.,  its  head  stops  changing).  The  steady-state  water  budget  in  this 
case  is  very  simple — all  supplies  and  demands  are  zero. 


The  above  example  illustrates  several  points  which  are  of  genera)  interest.  It 
is  apparent  that  an  initial  value  of  head  must  be  specified  throughout  the  solution 
region  if  the  simulation  problem  is  dynamic  (i.e.,  if  the  head  changes  over  time). 
The  time-varying  head  will  approach  a  steady-state  if  all  the  aquifer  inputs  are 
approximately  constant,  as  they  often  are  in  relatively  undeveloped  aquifers.  The 
manner  in  which  the  head  converges  to  steady-state  depends  strongly  on  the  initial 
head  distribution.  The  rate  at  which  this  convergence  takes  place  is  inversely 
proportional  to  the  storage  coefficient  (large  coefficients  give  slower  response  than 
small  coefficients). 


When  the  head  distribution  varies  with  location,  the  solution  of  the  groundwater 
flow  equation  becomes  more  complicated.  In  this  case,  the  auxiliary  constraints 
needed  to  insure  uniqueness  arc  specified  (for  all  times)  along  the  boundary  of  the 
solution  region.  Such  boundary  conditions  generally  take  the  following  forms: 

1.  SpeciGcd-Ilead  Boundary  Conditions  In  this  case,  the  head  h(x,  y,  t)  is  set  equal 
to  a  known  value  /q,(x,  y,  t)  over  some  specified  portion  of  the  boundary. 

2.  Specified-Flux  Boundary  Conditions  In  this  case,  the  normal  component  of  the 
vector  velocity  crossing  the  boundary  is  specified.  This  is  usually  accomplished 
by  imposing  the  constraint: 

-K  — =  qb(x,  y,  t )  (3-22) 

where  n  represents  the  distance  along  an  outward  pointing  vector  normal  to 
the  boundary  at  the  location  (x,y)  and  qb{x,y,t)  is  the  normal  component  of 
the  flux  into  the  aquifer  (in  units  l/t). 

Either  but  not  both  of  the  above  conditions  must  be  imposed  at  every  point 
along  the  boundary.  In  practice,  the  specified  flux  boundary  condition  is  usually 
only  used  when  the  boundary  is  known  to  be  a  barrier  to  flow  (i.e.,  when  qb{x,y,t) 
is  zero).  This  is  because  non-zero  boundary  flux  values  are  difficult  to  estimate 
from  available  data,  particularly  when  boundaries  are  located  far  from  monitored 
wells. 

The  second  example  in  this  section  considers  the  use  of  combined  head  and 
flux  boundary  conditions  in  a  hypothetical  problem  similar  to  the  one  examined  in 
Example  1. 


Example:  Boundary  Conditions 

Suppose  that  the  aquifer  from  the  preceding  example  is  not  totally  isolated 
but  is,  rather,  influenced  by  recharge  from  a  stream  that  runs  above  the  left-hand 
boundary.  This  effect  may  be  simulated  by  imposing  a  constant  head  boundary 
condition  (h(x,  y,  t)  =  ho)  along  the  entire  boundary.  Suppose,  in  addition,  that 
leakage  and  pumpage  are  negligible  but  that  the  infiltration  rate  has  a  uniform 
non-zero  value  of  q,o.  As  before,  the  transmissivity  and  initial  heads  are  assumed 
to  be  equal  to  T  and  ha  throughout  the  aquifer.  This  revised  problem  is  illustrated 
in  Figure  3-5. 

As  time  progresses,  the  infiltrating  water  from  above  will  begin  to  accumulate 
and  the  head  will  rise  above  its  uniform  initial  value  (except  along  the  fixed  head 
boundary).  Eventually,  the  internal  head  gradient  will  be  sufficiently  large  to  drive 
ail  the  infiltrating  water  out  across  the  left  side  of  the  aquifer.  The  head  will  then 
stop  changing  and  the  aquifer  will  be  in  steady-state.  The  steady-state  head  profile 
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FIGURE  3-5 

PLAN  AND  ELEVATION  OF  THE  AQUIFER  ANALYZED  IN  EXAMPLE  2 


may  be  found  by  setting  the  temporal  derivative  in  Equation  (3-14)  equal  to  zero, 
The  resulting  partial  differential  equation  is: 


d  „dh  d  dh 
dxTdx  +  dyTdy  +  <7‘°~ 


(3-23) 


Here  h  is  a  function  of  location  (x  and  y)  but  not  of  time. 

The  integration  of  Equation  (3-23)  may  be  simplified  by  noting  that  the  flow 
problem  defined  in  Figure  3-5  is  completely  uniform  in  the  y  direction,  i.e.,  there  is 
no  reason  why  the  head  should  vary  in  this  direction.  Consequently,  it  is  reasonable 
to  (tentatively)  assume  that  the  head  depends  only  on  the  x  coordinate.  This 
assumption  may  later  be  checked  by  substituting  the  resulting  head  distribution 
into  Equation  (3-23)  to  insure  that  it  is,  in  fact,  a  solution. 

If  h  depends  only  on  x,  dh/dy  is  zero  and  Equation  (3-14)  becomes  an  ordinary 
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differential  equation: 


££**-•  <3-24) 

This  may  be  integrated  once  to  give: 

dh 

Tdx=~q,oX  +  C'  (3~25) 

Here  C\  is  an  unknown  constant  of  integration.  Equation  (3-25)  may  be  solved  for 
dh/ dx  and  integrated  a  second  time  to  give  a  general  expression  for  the  head: 

y)  =  “  x2  +  ~~  +  C‘2  (3-26) 

where  Ci  is  a  second  integration  constant. 

Equation  (3-26)  gives  a  head  distribution  which  satisGcs  the  basic  groundwater 
equation  for  any  values  of  C\  and  C‘>.  This  non-uniqueness  is  eliminated  by 
imposing  the  boundary  conditions  defined  earlier.  If  there  is  no  tlow  across  the 
shaded  boundaries,  the  specific  discharge  across  these  boundaries  must  equal  zero. 
But  Darcy’s  law  (Eq.  3-6)  states  that  the  specific  discharge  along  the  right-hand 
boundary  is  given  by: 

,(L,y)  =  -K?^  (3-27) 

Since  the  aquifer  has  constant  depth  the  no-flow  condition  may  be  written  as: 


Tah(L,y) 

i  dx 

where  d  is  the  depth.  Equation  (3-25)  may  be  substituted  into  Equation  (3-28)  and 
solved  for  Cj  (note  that  i  is  equal  to  L  since  the  no-flow  condition  applies  at  the 
right-hand  end  of  the  aquifer): 

Cl  =  qi0L  (3-29) 

The  head  solution  now  has  only  one  arbitrary  constant  which  may  be  found  by 
noting  that  the  head  along  the  left-hand  boundary  (where  x  =  0)  is  equal  to  ho. 
Setting  x  —  0  in  Equation  (3-21)  gives: 

h{o,  y)  =  C2  =  h0  (3-30) 

The  unique  solution  obtained  subject  to  the  assumption  that  dh/dy  =  0  is  there¬ 
fore: 

h(x,y)  =  ~^x2+^~+h0  (3-31) 

Direct  substitution  reveals  that  this  expr-v  ,'jii  satisfies  Equation  (3-23).  It  clearly 
satisfies  the  left-  and  right-hand  boundary  conditions  since  it  was  constructed 
from  them.  It  satisfies  the  upper  and  lower  boundary  conditions  since  the  specific 
discharges  and  y  derivatives  across  these  boundaries  are  both  zero. 

The  head  and  velocity  distributions  for  this  example  are  plotted  vs.  x  in  Figure 
3-6.  Note  that  the  head  is  flat  at  the  right  end  where  the  flux  is  required  to  be  zero. 
The  head  at  the  left  end  is  fixed  at  the  value  ho,  as  specified  by  the  left-hand 
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boundary  condition.  The  velocity  vectors  point  in  the  x  direction  and  increase  in 
magnitude  from  right  to  left. 


This  example  clearly  shows  that  boundary  conditions  must  be  imposed  to 
give  a  unique  solution  to  the  steady-state  flow  equation.  If  a  dynamic  solution 
is  desired,  both  initial  conditions  and  boundary  conditions  must  be  imposed.  It 
should  be  noted  that  boundary  conditions  cannot  be  chosen  arbitrarily  -  they  must 
be  consistent  with  assumptions  made  in  the  flow  equation.  If,  for  example,  the 
aquifer  of  Figure  3-5  had  no-flow  boundaries  on  all  four  sides,  there  would  be  no 
way  for  the  inGltrating  water  to  leave.  This  water  would  accumulate  indefinitely, 
causing  the  head  to  rise  steadily  throughout  the  aquifer.  In  this  case,  there  is 
no  steady-state  solution,  i.c.,  there  is  no  head  distribution  which  can  satisfy  both 
Equation  (3-23)  and  the  imposed  boundary  conditions. 

Since  auxiliary  conditions  must  ultimately  be  estimated  from  limited  quantities 
of  field  data  the  model’s  boundaries  and  simulation  period  should  be  selected 
to  simplify  the  estimation  process  as  much  as  possible.  The  following  general 
guidelines  are  often  helpful: 

1.  Boundaries  should  be  laid  out,  whenever  possible,  along  flow  barriers  or  lines 
of  symmetry. 

2.  Specified  head  boundaries  should  pass  through  regions  where  the  head  is  rela¬ 
tively  constant  through  time.  It  is  best  if  these  boundaries  lie  near  monitored 
wells. 

3.  In  impact  studies  such  as  the  San  Andres-Glorieta,  where  the  effects  of  pumping 
are  of  primary  concern,  it  is  wise  to  locate  the  model’s  outer  boundaries  well 
beyond  the  region  likely  to  be  affected  by  the  pumpage.  Internal  boundaries 
maybe  used  to  describe  faults  (no- flow  boundaries)  or  surface  water  bodies  such 
as  lakes  or  perennial  streams  (constant  head  boundaries). 

4.  If  possible,  the  simulation  should  be  started  at  a  time  when  the  aquifer  is 
known  to  be  at  steady-state.  In  this  case,  a  steady-state  head  solution  may  be 
used  to  initialize  subsequent  dynamic  simulations. 

5.  If  a  steady-state  initialization  is  not  feasible,  the  simulation  should  be  started 
at  a  time  when  a  reasonable  number  of  reliable  well  observations  are  available 
for  defining  an  initial  condition. 

If  the  above  suggestions  are  followed,  the  model’s  data  requirements  can  be  reduced 
significantly  (see  Section  3.3.3). 

Most  groundwater  flow  problems  cannot  be  solved  using  the  analytical  in¬ 
tegration  procedures  outlined  in  the  examples.  These  procedures  require  assump¬ 
tions  of  uniformity,  symmetry,  and  geometric  regularity  which  rarely  hold  in  real 
aquifers.  Fortunately,  there  are  many  numerical  methods  for  integrating  the  flow 
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equations  when  boundary  conditions  and  aquifer  geometry  are  more  complex.  These 
methods  involve  repetitive  calculations  that  are  ideally  suited  for  electronic  com¬ 
puters.  Section  3.3  briefly  discusses  some  of  the  numerical  methods  and  computer 
programs  available  for  solving  realistic  groundwater  problems. 

3.2.4  Impact  Evaluation,  Drawdown,  and  Superposition 

The  procedure  for  simulating  a  vertically  averaged  head  distribution  for  the 
San  Andres- Glorieta  case  study  should  now  be  apparent.  First,  a  solution  region  is 
laid  out  and  boundary  conditions  are  defined.  If  the  simulation  is  dynamic,  an  initial 
head  distribution  must  also  be  specified.  Then  the  various  coefficients  and  source 
terms  in  the  flow  equation  are  estimated.  These  include  the  aquifer  parameters 
(transmissivity,  storage  coefficient,  and  leakance),  infiltration  and  pumpage  rates, 


and  adjoining  aquifer  heads.  Finally,  the  flow  equation  is  solved  for  the  hydraulic 
head  distribution.  The  hydrologic  impacts  of  different  management  strategies  are 
evaluated  by  comparing  the  head  solutions  obtained  with  appropriate  pumping 
rates. 

Hydrologic  impact  evaluation  is  concerned  primarily  with  the  change  in  ground- 
water  elevation  which  results  from  a  change  in  some  spccifed  (or  “nominal”)  pump¬ 
ing  strategy.  The  nominal  strategy  may  be  selected  in  many  ways  -  it  could  repre¬ 
sent  undeveloped  conditions  (no  pumpage),  it  could  be  a  continuation  of  current 
pumping  levels,  or  it  could  be  a  “low  growth”  option  based  on  a  slowly  increasing 
pumpage  rate.  The  head  associated  with  the  nominal  pumping  strategy  may  be 
called  the  nominal  head  while  the  head  associated  with  some  modified  pumping 
strategy  (usually  more  pumpage)  may  be  called  the  modified  head.  It  is  customary 
in  impact  studies  to  use  the  nominal  head  as  a  reference  for  computing  drawdown, 
which  is  then  defined  as: 

d(x,  y,  t )  =  hm(x,  y,  t)  -  hn( x,  y,  t)  (3-32) 

Here  the  subscripts  n  and  m  represent  nominal  and  modified  heads,  respectively. 

These  concepts  are  illustrated  in  Figure  3-7,  which  shows  nominal  and  modified 
heads  and  drawdown  for  a  simple  impact  evaluation  problem.  The  nominal  pump¬ 
ing  strategy  consists  of  pumpage  from  a  single  well.  This  pumpage  results  in  a 
general  decline  in  water  levels,  as  indicated  in  Figure  3-7a.  The  modified  strategy 
illustrated  in  Figure  3-7b  includes  a  second  well  which  further  accelerates  the  rate 
of  decline.  The  drawdown  curves  plotted  in  Figure  3-7c  are  obtained  by  differencing 
nominal  and  modified  heads. 

The  situation  pictured  in  Figure  3-7  may  be  described  mathematically  with 
two  vertically  averaged  flow  equations,  one  for  each  pumping  strategy: 


Nominal  pumping  strategy 

odhn  _  d  dhn 
dt  dx  dx 


dy  dy 


—  9.  +  Qn 


+  Lu(hv  —  hn )  +  Li(h[  —  /in) 


(3-33) 


h»{x,y,o)  =  h0(x,y) 


initial  conditon 


hn{x,y,t)  =  hb{x,y,t) 
-K^~{x,y,t)  =  qb(x,y,t) 


head  boundary  condition 


flux  boundary  condition 
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(3-34) 


Modified  pumping  strategy 
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head  boundary  condition 
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The  terms  preceded  by  A  represent  the  changes  in  pumpage  and  recharge  associated 
with  the  modified  pumping  schedule  (recharge  might  change  if  some  of  the  addi¬ 
tional  water  pumped  eventually  infiltrates  back  to  the  aquifer).  Note  that  the 
auxiliary  conditions  and  the  heads  in  adjoining  aquifers  ( hu  and  hi)  are  assumed  to 
be  the  same  for  the  nominal  and  modified  pumpage  strategies.  This  implies  that 
both  strategies  start  with  the  same  initial  head  distribution.  It  also  implies  that 
the  boundary  conditions  are  unaffected  by  pumpage. 

If  the  aquifer  is  confined  or  if  it  is  unconfincd  but  head  variations  are  small 
compared  to  the  saturated  flow  depth,  the  aquifer  parameters  S,  T,  Lu,  and  L\ 
have  the  same  (constant)  values  in  each  equation.  In  this  case,  the  problem  is  linear 
and  the  principle  of  superposition  holds  (Freeze  and  Cherry,  1979).  The  modified 
pumpage  equation  may  then  be  subtracted  from  the  nominal  pumpage  equation  to 
give  the  following  expression  for  drawdown: 


Drawdown 


dd  d  _<9d  d  dd 

S  -r-  -  j-r - - —T  —  =  A +  A Q  -  Lud  -  Lid 

at  ox  ox  oy  ay 


(3-35) 


d(x,  y,  o ) 
d{x,  y,  t ) 
dd, 

-zz(x>y>  *) 


initial  condition 

drawdown  boundary  condition 

gradient  boundary  condition 


This  important  equation  reveals  that  drawdown  depends  only  on  the  aquifer  para¬ 
meters  and  the  assumed  changes  in  pumpage  and  recharge.  The  initial  and  bound¬ 
ary  conditions  for  the  drawdown  equation  arc  identically  zero  (and  so  are  known 
perfectly)  and  the  nominal  pumpage  and  infiltration  drop  out. 

The  result  derived  above  suggests  that  the  inputs  to  a  vertically  averaged  ground- 
water  model  can  be  divided  into  two  different  categories  --  primary  inputs  which 
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affect  hot h  head  and  drawdown  predictions  and  secondary  inputs  which  affect  only 
head.  Those  categories  are  defined  as  follows: 

Pri m ary  Inputs 

a)  Time-invariant  Aquifer  Parameters 

Transmissivity 
Storage  coefficient 
Leakancc 

b)  Postulated  Changes  in  Groundwater  Fluxes 

Purnpagc  change 
Recharge  change 

Secondary  Inputs 

a)  Auxiliary  Conditions 

Initial  head 
Boundary  head  or  flux 
(Adjoining  aquifer  heads) 

b)  Nominal  Values  of  Groundwater  Fluxes 

Nominal  pumpage  and  discharge 
Nominal  infiltration 

Since  drawdown  is  of  particular  interest  in  impact  evaluations,  the  primary  inputs 
deserve  the  most  time  and  attention  during  the  modeling  process.  Secondary  inputs 
are  important  in  several  respects  mentioned  later  but  clearly  do  not  have  as  much 
influence  on  the  overall  conclusions  of  the  impact  analysis. 

It  might  seem  reasonable  to  forget  head  simulation  altogether  and  simply  solve 
the  drawdown  equation  directly.  This  is  a  legitimate  approach  which  is  frequently 
used  (see  Chapter  4  for  an  example).  Drawdown  simulation  does,  however,  have 
some  practical  limitations  which  are  worth  noting. 

It  is  important  to  remember  that  the  drawdown  derivation  is  based  on  the 
assumption  of  linearity.  Strictly  speaking,  the  groundwater  equation  is  linear  only 
if  the  aquifer  is  completely  confined.  If  unconfined  effects  are  significant,  the 
transmissivity  and  storage  coefficient  both  depend  on  head  and  the  superposition 
operation  needed  to  derive  the  Equation  (3-35)  cannot  be  carried  out.  This  implies 
that  drawdown  models  should  not  be  used  in  areas  where  significant  depletion  or 
dewatering  may  take  place. 

Although  drawdown  is  the  primary  measure  of  the  hydrologic  impact  of  de¬ 
velopment  it  docs  not  tell  the  whole  story.  In  some  situations,  it  is  important  to 


know  whether  the  groundwater  elevation  will  fall  below  a  specified  level,  such  as  the 
top  of  the  aquifer  or  the  inlet  of  an  existing  well.  Such  questions  can  be  answered 
only  with  a  head  simulation  which  is  referred  to  a  fixed  datum.  The  secondary 
inputs  required  to  simulate  head  but  neglected  in  drawdown  simulation  provide  the 
information  needed  to  locate  groundwater  levels  relative  to  such  a  datum. 

There  arc  clearly  some  advantages  to  the  straightforward  approach  of  simulat¬ 
ing  both  nominal  and  modified  heads.  The  drawdown  predictions  obtained  by 
differencing  these  heads  will  be  the  same  as  those  obtained  by  solving  the  draw¬ 
down  equation  directly  and  additional  information  will  be  gained  on  absolute  water 
levels.  Comparisons  of  model  predictions  and  field  observations  arc  often  easier  if 
the  head  simulation  approach  is  used.  The  only  real  disadvantage  of  this  approach 
is  the  extra  effort  required  to  estimate  secondary  model  inputs.  If  information  on 
these  inputs  is  very  limited  a  simulation  based  on  the  drawdown  equation  may  be 
the  best  alternative. 


3.3  MODEL  APPLICATION 


3.3.1  Selection  of  a  Computer  Program 

Nearly  all  computer  programs  used  to  simulate  the  hydrologic  impacts  of 
groundwater  development  arc  based  on  the  same  general  concepts  -the  principles  of 
mass  continuity  and  Darcian  flow  discussed  in  Section  3.2.  The  primary  differences 
between  available  programs  are  the  dimensionality  of  the  governing  flow  equa¬ 
tion  (either  the  full  three-dimensional  equation  or  the  two-dimensional  vertically 
averaged  equation  is  usually  used)  and  the  numerical  procedure  used  to  integrate 
this  equation. 

Most  groundwater  modeling  studies  rely  on  two-dimensional  (vertically  averaged) 
computer  models,  although  three-dimensional  models  arc  being  used  with  increasing 
frequency.  Three-dimensional  analyses  are  generally  required  if  groundwater  is  be¬ 
ing  pumped  from  several  aquifers  which  are  hydrologically  connected.  In  such  cases, 
ihe  simple  leakage  relationship  used  in  the  two-dimensional  vertically  averaged  flow 
equation  may  not  be  adequate.  Although  three-dimensional  models  can,  in  prin¬ 
ciple,  provide  a  more  realistic  description  of  complex  flow  patterns  they  require 
more  input  data  and  are  more  expensive  to  run  than  their  two-dimensional  coun¬ 
terparts.  Some  three-  dimensional  inputs,  such  as  the  vertical  component  of  the 
hydraulic  conductivity,  are  very  difficult  to  estimate  from  Geld  observations.  The 
actual  bcneGts  of  using  a  three-dimensional  model  depend  on  the  availability  of 
reliable  data  as  well  as  on  the  complexity  of  the  groundwater  system. 

Several  different  numerical  procedures  are  available  for  solving  either  the  two 
or  three-dimensional  groundwater  flow  equation.  The  most  popular  arc  the  Gnite 
difference  and  Gnite  element  techniques  described  in  Section  3.3.2.  Either  of  these 
will  give  acceptable  results  for  most  problems.  The  Gnite  element  method  is 
somewhat  more  convenient  to  use  with  solution  regions  which  have  irregularly- 
shaped  boundaries.  The  Gnite  difference  technique  is  often  less  expensive  and  is 
easier  to  understand  and  explain. 

Although  dimensionality  and  solution  technique  are  important,  there  are  also 
a  number  of  other  factors  to  be  considered  when  selecting  a  groundwater  program. 
The  program  should  be  well-documented,  preferably  with  a  user’s  manual  which 
includes  several  solved  examples.  Data  entry  should  be  convenient  and  the  program 
should  be  “portable”  between  different  computers.  The  program’s  outputs  should 
be  easy  to  understand  and  should  include  high-resolution  contour  plots  which  are 
able  to  display  both  measured  and  simulated  heads.  Finally,  the  program  should 
have  a  proven  record  of  successful  application  in  “real  world”  situations. 

Perhaps  the  most  widely  used  groundwater  flow  models  available  to  the  general 
public  are  the  U.S.  Geological  Survey  (USGS)  finite  difference  models  described  in 
Trescott,  Pindcr,  and  Larson  (1976)  and  in  Trcscott  (1975).  Several  state  agencies 
such  as  the  California  Department  of  Water  Resources  and  the  Kansas  Geological 
Survey  have  also  developed  models  which  arc  in  the  public  domain.  The  U.S.  Army 
Corps  of  Engineers’  Hydrologic  Engineering  Center  has  assembled  a  package  of 


model  application  programs  designed  to  be  used  with  the  USGS  two-dimensional 
finite  difference  model  (MFC,  1983).  This  package  includes  a  program  for  construct¬ 
ing  model  input  files  and  a  program  for  plotting  model  predictions.  A  somewhat 
out-of-date  but  useful  review  of  publicly  available  groundwater  programs  is  given  in 
the  American  Geophysical  Union’s  monograph  on  groundwater  modeling  (Bachmat 
et  al.,  1980). 

Although  there  are  many  groundwater  modeling  programs  to  choose  from, 
the  first-time,  or  infrequent  modeler  would  probably  do  best  to  stay  with  the  well- 
documented  and  thoroughly-tested  models  of  the  USGS.  This  choice  should  allow 
the  user  to  concentrate  on  problem  formulation  and  data  analysis  rather  than  the 
mechanics  of  computer  programming. 

3.3.2  Spatial  and  Temporal  Discretization 

Most  numerical  techniques  for  integrating  the  groundwater  flow  equation  con¬ 
vert  the  original  partial  differential  equation  into  a  large  set  of  algebraic  equations 
which  may  be  readily  solved  on  a  digital  computer.  This  conversion  process  requires 
the  simulation  problem  to  be  discretized  in  both  space  and  time.  That  is,  the  solu¬ 
tion  region  is  divided  into  a  number  of  discrete  subregions  (or  elements)  and  the 
simulation  period  is  divided  into  a  number  of  discrete  time  intervals.  The  partial 
differential  equation  is  then  used  to  derive  one  or  more  algebraic  equations  for  each 
element  and  each  time  interval.  The  unknowns  in  these  algebraic  equations  are 
“average”  heads  which  approximate  the  exact  solution  at  a  given  time  and  location. 
As  the  discretization  is  made  more  detailed,  the  approximate  solution  converges  to 
the  exact  solution  everywhere. 

There  are  a  number  of  methods  for  discretizing  the  solution  regions  for  ground- 
water  problems.  Although  these  methods  appear  to  be  quite  different  they  share 
a  number  of  basic  concepts.  Figure  3-8a  shows  a  typical  irregularly-shaped,  two- 
dimensional  solution  region.  This  could  be  a  vertically  averaged  representation  of 
a  complete  aquifer  or  of  some  portion  of  an  aquifer.  Two  alternative  methods  of 
dividing  up  or  discretizing  the  solution  region  are  illustrated  in  Figures  3-8b  and 
3- 8c. 

Figure  3-8b  superimposes  a  grid  constructed  of  unevenly  spaced  parallel  lines 
on  top  of  the  solution  region.  The  intersections  of  these  lines  are  called  nodes 
and  the  rectangular  region  surrounding  each  node  is  called  a  grid  cell.  Most 
computer  models  which  use  this  type  of  discretization  scheme  assume  that  the  head 
everywhere  within  each  grid  cell  is  equal  to  the  head  at  the  corresponding  node 
(other  assumptions  are  possible).  If  a  rectangular  spatial  discretization  technique 
is  used,  the  flow  equation  may  be  converted  into  a  set  of  N  ordinary  differential 
equations,  where  N  is  the  number  of  interior  nodes  in  the  grid.  The  heads  or  head 
gradients  at  the  boundary  nodes  are  obtained  from  boundary  conditions.  Each 
ordinary  differential  equation  depends  only  on  time  and  may  be  further  discretized 
using  methods  discussed  below. 


B)  Actual  aquifen  outline 
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b)  Rectangular  simulation  qrid 
(e-9-  finite  difference) 
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FIGURE  3-8 

SOME  TYPICAL  SPATIAL  DISCRETIZATION  SCHEMES 


Figure  3- 8c  shows  an  irregular  discretization  network  which  divides  the  solution 
region  into  elements  with  curved  sides.  In  this  case,  the  head  within  each  element 
is  assumed  to  be  a  weighted  sum  of  the  heads  at  the  surrounding  node  points.  As 
before,  the  How  equation  is  converted  into  a  set  of  /V  ordinary  differential  equations, 
where  N  is  the  number  of  interior  nodes.  The  rectangular  grid  approach  is  typically 
used  with  finite  difference  models  such  as  the  USCS  finite  difference  model  selected 
for  the  San  Andres-Clorieta  study.  The  irregular  curve-sided  approach  is  typically 
used  with  finite  element  models. 

The  ordinary  differential  equations  generated  by  spatial  discretization  are  gen¬ 
erally  discret  ized  in  time  using  a  method  similar  to  the  one  illustrated  in  Figure  3-9. 
The  input  variables  and  head  derivatives  appearing  in  the  equation  are  assumed 
to  be  constant  over  each  time  interval  and  to  jump  abruptly  at  the  end  of  the 
interval.  Since  the  head  derivative  is  assumed  constant,  the  head  itself  is  assumed 
to  vary  linearly  over  the  time  interval.  When  temporal  discretization  is  applied,  the 
each  ordinary  differential  equation  is  replaced  by  a  set  of  algebraic  equations  whose 
unknowns  are  the  nodal  heads  at  the  end  of  each  time  interval.  Note  that  the  head 
at  the  beginning  of  the  first  time  interval  is  obtained  from  the  initial  condition. 

The  process  of  laying  out  a  discretized  grid  or  network  for  a  particular  model¬ 
ing  problem  depends  both  on  the  type  of  mode!  selected  and  on  the  specific  require¬ 
ments  of  the  problem.  Program  user  manuals  usually  give  guidelines  for  network 
construction.  Generally  speaking,  the  construction  of  finite  difference  grids  is  more 
mechanical  and  straightforward  than  the  construction  of  finite  element  networks. 
The  smallest  cells  in  a  finite  difference  grid  are  usually  located  in  the  area  where 
impacts  are  of  most  concern.  In  most  cases  this  is  the  region  surrounding  a  well 
field.  Grid  spacing  usually  increases  in  both  directions  away  from  this  region. 

Finite  element  networks  are  usually  less  detailed  than  finite  difference  grids 
for  the  same  aquifer.  The  clement  sides  in  these  networks  typically  coincide  with 
geological,  hydrological,  or  institutional  boundaries  which  are  relevant  to  the  model¬ 
ing  study.  A  given  element  might,  for  example,  be  defined  as  the  area  within  a 
particular  irrigation  district  which  is  devoted  to  grazing  land  and  characterized  by 
sandy  soils.  Such  areas  can  be  identified  by  overlaying  transparent  maps  showing 
land  uses,  soil  types,  administrative  regions  and  other  relevant  spatial  features. 
The  superimposed  composite  of  all  of  these  maps  usually  gives  a  good  preliminary 
network.  It  is  easy  to  either  refine  or  simplify  this  network  if  additional  information 
becomes  available  later.  This  flexibility  is  one  of  the  most  attractive  features  of  the 
finite  element  approach. 

3.3.3  Input  Estimation 

Input  estimation  is  probably  the  most  important  and  most  neglected  single  task 
in  the  modeling  process,  particularly  in  groundwater  applications  where  only  a  few 
of  the  relevant  variables  are  directly  observable.  This  is  clearly  illustrated  in  the 
detailed  analysis  of  the  San  Andrcs-Glorieta  modeling  studies  presented  in  Chapter 
4.  The  need  to  estimate  inputs  throughout  the  solution  region  and  over  a  simulation 
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FIGURE  3-9 

SOME  TYPICAL  TEMPORAL  DISCRETIZATION  SCHEMES 


period  extending  many  years  into  the  future  forces  the  modeler  to  extrapolate  and 
generalize  from  the  limited  data  available.  This  inevitably  introduces  subjectivity 
and  uncertainty  into  the  modeling  process. 

It  is  convenient  to  divide  this  discussion  of  input  estimation  techniques  into 
two  sections—  one  dealing  with  aquifer  parameters  (primary  inputs)  and  one  dealing 
with  auxiliary  conditions  and  nominal  pumpage  and  recharge  (secondary  inputs). 
Although  the  emphasis  here  is  on  vertically  averaged  models  of  hydraulic  head,  most 
of  the  techniques  considered  are  also  applicable  to  drawdown  models,  provided  the 
superposition  assumption  holds  (see  Section  3.2.4).  Some  practical  applications  of 
these  techniques  arc  discussed  in  Chapter  4. 


i)  aquifer  parameters 


Simulations  of  hydraulic  head  and  drawdown  require  estimates  of  transmis¬ 
sivity,  storage  coefficient,  and  leakancc  at  each  node  or  element  in  the  model  net¬ 
work  (input  values  need  not,  of  course,  be  different  at  every  node  or  clement). 
These  physical  parameters  generally  vary  throughout  the  aquifer  variations  tend 
to  be  large  in  regions  which  are  geologically  complex  and  small  in  regions  which 
are  homogeneous.  This  is  illustrated  in  a  simple  way  in  Figure  ?  10b,  which  shows 
a  plot  of  transmissivity  vs.  distance  along  a  transect  in  a  hypothetical  aquifer. 
Although  the  regional  average  is  approximately  1000  ff  /day,  local  values  based, 
for  example,  on  soil  samples  may  vary  from  100  ft“/day  to  as  high  as  10,000  fl2/day. 
This  type  of  variability  is  typical  of  the  data  available  for  the  San  Andrcs-Glorieta. 

The  aquifer  parameter  estimates  produced  by  most  estimation  techniques  are 
averages  which  apply  over  a  characteristic  range  or  length  scale.  This  is  illustrated 
in  Figure  3- 10a  where  the  ranges  for  estimates  derived  from  soil  samples,  pump 
tests,  and  regional  model  calibration  are  compared.  The  estimates  produced  by 
each  of  these  techniques  are  indicated  by  appropriate  horizontal  lines  in  Figure 
3- 10b. 

Soil  samples  and  well  logs  obtained  during  well  drilling  give  a  very  localized  in¬ 
dication  of  the  grain  size  distribution  in  a  particular  portion  of  an  aquifer.  Various 
empirical  formulas  are  available  for  estimating  hydraulic  conductivity  and  trans¬ 
missivity  from  this  distribution  (sec  Freeze  and  Cherry,  1979).  These  should  be 
regarded  with  skepticism  and  used  only  in  the  absence  of  any  better  alternatives. 

Pump  tests  are  probably  the  most  popular  method  for  estimating  aquifer 
transmissivities  and  storage  coefficients,  partly  because  they  give  a  good  indication 
of  a  well’s  producing  capacity.  The  pump  test  approach  is  an  example  of  an  inverse 
method  of  parameter  estimation.  The  selected  well  is  pumped  for  a  specified 
period  and  water  levels  are  observed  during  both  the  drawdown  and  recovery 
process.  A  simplified  (radially  symmetric)  groundwater  equation  is  then  solved 
for  the  transmissivity  and  storage  coefficient  which  best  reproduce  the  water  level 
history  observed  during  the  test.  This  process  has  been  standardized  by  Theis 
(1935),  Jacob  (1940),  Hantush  (1960)  and  others  (see  Freeze  and  Cherry,  1979,  for 
a  summary). 

The  range  over  which  pump  test  results  apply  depends  on  a  number  of  factors 
including  the  duration  of  the  test  -  brief  tests  only  produce  drawdowns  close  to 
the  pumped  well  while  longer  tests  affect  a  wider  area.  For  aquifers  as  large  as  the 
San  Andrcs-Glorieta,  pump  test  estimates  should  probably  be  interpreted  as  point 
observations  (although  they  are  admittedly  less  localized  than  soil  samples). 
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FIGURE  3-10 

SPATIAL  SCALES  OF  DIFFERENT  PARAMETER  ESTIMATION  TECHNIQUES 
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FIGURE  3-11 

ALTERNATIVE  WAYS  TO  EXTRAPOLATE  LOCALIZED  WELL  OBSERVATIONS 


If  pump  test  results  are  to  be  used  to  define  aquifer  parameters  throughout 
the  solution  region  (i.e.,  at  all  nodes  or  elements)  they  must  be  extrapolated. 
Two  possible  extrapolation  methods  are  illustrated  in  Figure  3-11.  Figure  3- 11a 
shows  a  hypothetical  aquifer  which  contains  five  observation  wells.  Transmissivity 
estimates  obtained  from  pump  tests  at  these  wells  are  indicated  next  to  each  well 
symbol.  Figure  3-llb  shows  a  blocked  extrapolation  technique  which  allocates 
each  transmissivity  value  to  a  large  geologically  homogeneous  area  surrounding  a 
particular  well  (or  wells).  Figure  3-1  lc  shows  an  extrapolation  technique  based 
on  a  contour  plot  of  the  well  values.  The  contours  may  be  interpolated  to  give  a 
continuously-varying  distribution  of  transmissivity.  This  approach  is  particularly 
appropriate  in  geologically  heterogeneous  aquifers  such  as  the  San  Andres-Glorieta. 

Regional  model  calibration  is  an  estimation  procedure  which  avoids  some  of 
the  disadvantages  associated  with  extrapolation.  This  procedure  is  an  inverse 


technique  similar  in  concept  to  pump  test  analysis  but  quite  different  in  its  practical 
application.  Regional  calibration  uses  the  complete  vertically  averaged  flow  model 
to  simulate  heads  throughout  the  aquifer  over  a  historical  observation  period. 
The  unknown  aquifer  parameters  are  adjusted,  either  by  trial-and-crror  or  with 
an  optimization  algorithm,  until  a  “best  fit”  is  obtained  between  the  historically 
observed  heads  and  the  heads  simulated  by  the  model.  This  parameter  adjustment 
process  is  the  counterpart  to  the  type-curve  fitting  procedures  used  to  obtain 
estimates  from  pump  test  measurements. 

It  is  important  to  note  that  the  storage  coefficient  can  be  estimated  with  a 
regional  calibration  approach  only  if  heads  during  the  historical  period  are  varying 
(i.c.,  are  not  in  steady-state).  Otherwise,  the  temporal  head  derivative  is  zero  and 
the  storage  coefficient  has  no  effect  on  the  model’s  predictions.  In  practice,  it 
is  probably  best  to  accept  this  and  stay  with  a  steady-state  simulation  whenever 
possible.  This  is  because  initial  condition  errors  introduced  in  a  dynamic  simulation 
could  easily  outweigh  any  information  gained  about  the  storage  coefficient. 

The  major  advantage  of  regional  model  calibration  is  its  ability  to  provide 
estimates  of  the  spatially  discretized  (regionally  averaged)  parameters  used  in  the 
model.  These  estimates  do  not  have  to  be  extrapolated  or  generalized  in  any  way  - 
they  can  be  input  to  the  model  as  is.  The  major  disadvantage  of  regional  calibration 
is  its  dependence  on  the  accuracy  of  the  input  values  entered  for  historical  pumpage, 
recharge,  and  boundary  conditions.  If  these  inputs  are  incorrect,  the  estimated 
aquifer  parameters  may  be  quite  far  from  the  true  values. 

When  all  the  alternatives  are  considered,  the  best  approach  to  aquifer  parameter 
estimation  seems  to  be  a  combination  of  the  pump  test  and  regional  calibration  tech¬ 
niques.  This  can  be  accomplished  by  starting  a  regional  calibration  with  parameter 
estimates  obtained  from  pump  tests.  Subscqcnt  parameter  adjustments  should  be 
constrained  so  that  they  improve  upon  but  do  not  deviate  too  far  from  the  initial 
pump  test  values.  This  approach  seems  to  work  well  in  practice,  probably  because 
it  uses  all  available  sources  of  information  in  an  efficient  way. 


ii)  secondary  inputs 


The  secondary  inputs  required  for  a  simulation  of  hydraulic  head  are  initial 
heads  (if  the  simulation  is  dynamic),  boundary  conditions  (including  the  heads  in 
adjoining  aquifers),  and  nominal  pumpage  and  recharge,  'flic  estimation  problem 
is  clearly  much  easier  if  the  simulation  is  steady-state  since  boundary  conditions, 
pumpage,  and  recharge  arc  then  all  constants.  Fixed  boundary  head  or  fluxes 
are  usually  estimated  from  water  level  observations  collected  in  the  interior  of  the 
solution  region  where  pumpage  is  taking  place.  These  must  be  extrapolated  out  to 
the  boundaries  with  a  technique  similar  to  the  one  illustrated  in  Figure  3-llc. 

If  a  dynamic  simulation  is  required  it  is  best  to  start  computing  at  a  time 
when  the  aquifer  is  in  steady-state.  In  this  case,  the  initial  heads  can  be  derived 
directly  from  a  steady-state  solution.  Otherwise  these  heads  must  be  estimated  by 
extrapolating  well  observations  taken  at  the  initial  time.  This  can  be  a  significant 
source  of  prediction  error,  particularly  during  the  beginning  of  the  simulation 
period,  if  the  model  is  forced  to  adjust  to  a  physically  unrealistic  initial  head 
distribution. 

When  the  boundary  conditions  used  in  a  dynamic  head  simulation  are  time- 
varying  the  input  estimation  problem  becomes  very  difficult,  if  not  impossible.  This 
can  be  avoided  if  the  model’s  boundaries  arc  located  along  lines  where  the  head  or 
flux  are  reasonably  constant,  as  suggested  in  Section  3.2.3.  Flow  barriers  and  lines 
of  symmetry  should,  of  course,  be  used  to  deflne  boundaries  whenever  possible. 

One  might  expect  that  historical  pumping  rates  could  be  readily  estimated 
from  water  use  or  power  consumption  records  or  indirectly  derived  from  population 
and  land  use  data.  Although  these  information  sources  arc  all  helpful,  they  do 
not  necessarily  provide  a  completely  accurate  record  of  aquifer  pumping.  The 
estimation  problem  is  complicated  considerably  when  major  water  supply  wells 
arc  completed  in  more  than  one  aquifer  or  when  pumping  is  sporadic  (depending, 
for  example,  on  crop  demands  and  weather).  This  is,  for  example,  the  situation 
encountered  in  the  San  Andres-Glorieta.  It  is  often  useful  to  construct  a  set  of 
surface  water  budgets  for  major  land  uses  (particularly  in  agricultural  regions)  so 
that  pumpage  estimates  can  be  checked  for  consistency  with  other  surface  water 
data. 

Although  pumpage  can  often  be  estimated  for  a  single  well  or  a  small  well 
field,  it  should  be  remembered  that  the  spatial  resolution  of  a  discretized  model 
is  limited  by  the  size  of  its  elements  or  grid  cells.  Pumpage  from  a  point  located 
within  a  given  clement  is  effectively  spread  over  the  element.  If  the  modeler  wishes 
to  localize  pumpage  precisely  he  must  make  the  elements  in  the  vicinity  of  the  well 
very  small. 

Recharge  to  confined  aquifers  such  as  the  San  Andres-Glorieta  occurs  primarily 
near  outcrop  areas  where  the  aquifer  is  exposed  to  the  surface  or  is  overlain  by 
permeable  strata.  The  net  long-term  (c.g.,  annual)  recharge  may  be  estimated  by 
subtracting  evapotranspiration  and  runoff  from  total  precipitation  (including  snow- 


melt).  Regional  evapotranspiration  rates  may  he  estimated  from  various  empiri¬ 
cal  formulas  or  sometimes  derived  from  field  observations.  Average  precipitation 
data  arc  usually  readily  available  from  the  National  Weather  Service  or  from  local 
sources.  Runoff  may  be  available  from  stream  records  or,  alternatively,  may  be 
estimated  using  a  variety  of  empirical  techniques  (sec  Linsley  et  ai. ,  1977).  One  of 
the  greatest  sources  of  uncertainty  in  recharge  estimation  is  the  size  of  the  recharge 
area.  This  must  be  deduced  from  geological  observations  which  arc  often  not  very 
extensive  or  specific. 

Nominal  pumpage  and  recharge  for  the  future  can  be  projected  or  postulated 
in  a  number  of  ways.  Both  pumpage  and  recharge  could  be  assumed  to  remain  fixed 
at  current  levels.  Or  pumpage  could  be  held  at  current  levels  and  recharge  could 
be  varied  according  to  some  specified  climatic  sequence.  This  approach  allows  the 
modeler  to  investigate  the  impact  of  prolonged  droughts  or  wet  periods.  If  the  region 
being  studied  is  growing  rapidly  it  may  be  reasonable  to  assume  that  the  nominal 
pumpage  increases  over  time.  It  should  be  apparent  that  the  process  of  estimating 
nominal  pumpage  and  recharge  for  the  future  depends  on  assumptions  made  about 
water  supply  management  policies,  regional  growth,  and  climatic  conditions.  These 
assumptions  should,  of  course,  be  carefully  spelled  out  when  the  model’s  results  are 
reported. 

3.3.4  Accuracy  Evaluation 

It  should  now  be  apparent  that  it  is  difficult  to  estimate  the  inputs  of  a  “real- 
world”  groundwater  model  without  introducing  significant  error  at  one  point  or 
another.  The  following  list  identifies  some  of  the  critical  stages  where  errors  can 
arise  in  a  study  such  as  the  San  Andres- Glorieta: 

1.  The  entire  modeling  process  depends  on  the  assumptions  made  in  the  govern¬ 
ing  equations.  If  the  aquifer  is  assumed  to  be  confined  when  it  is  actually 
unconfined,  or  in  steady-state  when  its  heads  are  actually  changing,  estimation 
and  prediction  errors  will  result. 

2.  The  appropiate  values  for  model  inputs  depend  on  the  discretization  method 
selected.  Each  model  or  elemental  input  represents  an  average  over  an  extended 
region.  If  the  model  network  is  too  coarse,  these  averaged  inputs  may  not 
properly  represent  heterogeneities  in  the  actual  aquifer.  If  it  is  very  fine, 
computational  requirements  may  become  prohibitive. 

3.  Most  groundwater  inputs  must  be  estimated  from  a  small  number  of  measure¬ 
ments  taken  at  scattered  wells.  The  extrapolation  required  to  extend  these 
point  measurements  over  a  wide  area  can  introduce  significant  error. 

4.  Well  water  level  measurements  can  sometimes  be  misleading,  particularly  if  the 
well  is  multiply  completed  or  defective  in  some  way.  Water  level  measurement 
errors  affect  the  validity  of  pump  test  analyses,  extrapolation,  and  other  aspects 
of  the  estimation  process. 


5.  Traditional  pump  lost  analyses  assume  constant  pumpage,  homogeneous  aquifer 
properties,  negligible  well  loss,  complete  well  penetration,  and  regular  aquifer 
geometry.  Parameter  estimates  based  on  such  simplified  analyses  may  be  in¬ 
accurate  in  practical  applications.  Although  pump  test  estimation  errors  can 
be  reduced  if  more  sophisticated  (and  complex)  analytical  methods  are  used, 
they  can  never  be  completely  eliminated. 

6.  Aquifer  parameter  estimates  obtained  from  a  regional  model  calibration  are 
highly  dependent  on  the  values  used  for  inputs  such  as  historical  pumpage, 
recharge,  and  boundary  conditions.  If  these  inputs  arc  incorrect,  the  resulting 
parameter  estimates  may  be  inaccurate. 

In  a  realistic  groundwater  study  it  is  likely  that  one  or  more  of  these  errors  will  be 
significant.  It  is  important  to  have  some  idea  of  the  effects  these  errors  will  have  on 
model  predictions,  particularly  if  the  model’s  results  will  be  used  to  guide  policy 
decisions. 

A  view  of  input  estimation  and  subsequent  predictive  simulations  which  is 
particularly  relevant  to  the  San  Andres-Glorieta  application  is  illustrated  in  Figure 
3-12.  Figure  3-12a  shows  a  typical  time-varying  input  (e.g.,  recharge)  plotted 
over  the  past,  present,  and  future.  Measured  values  available  at  discrete  times 
in  the  past  are  used  to  estimate  the  actual  historical  recharge  (solid  line)  and  to 
project  recharge  trends  for  the  future.  Note  that  both  the  measured  values  and 
the  smoothed  estimates  differ  from  the  actual  recharge,  which  is  unknown  to  the 
modeler.  The  shaded  area  represents  the  input  estimation  error. 

Figure  3- 12b  shows  a  time- varying  observable  output  (e.g.,  head)  generated  by 
a  model  which  uses  the  estimated  input  record  plotted  in  Figure  3- 12a.  Here  again, 
observations  are  available  at  discrete  times  during  the  historical  period  and  the 
shaded  area  represents  the  model’s  prediction  error.  As  can  be  seen,  the  historical 
prediction  errors  do  not  necessarily  reveal  anything  about  the  model’s  prediction 
accuracy  in  the  future.  It  should  be  recognized  that  measurements  of  the  input 
and  output  variables  do  not  necessarily  coincide  with  the  unknown  actual  values. 
Accuracy  analyses  which  assume  that  all  measurements  are  perfect  are  unrealistic 
and  potentially  misleading. 

There  are  two  approaches  to  the  evaluation  of  model  accuracy,  one  which  relies 
on  “goodness  of  fit”  comparisons  between  model  predictions  and  field  observations 
and  one  which  relies  on  sensitivity  analysis.  The  objective  in  either  case  is  to 
obtain  a  quantitative  measure  of  the  model’s  prediction  error.  The  prediction  error 
is  unknown  by  definition  (otherwise  it  would  not  be  an  error)  but  it  can  still  be 
bounded  or  measured  in  a  probabilistic  sense.  Bounds  or  confidence  limits  on 
prediction  errors  provide  an  intuitive  and  helpful  guide  for  decision  making.  In  the 
San  Andres-Glorieta  case  study,  for  example,  “best”,  “worst”,  and  “most  likely” 
predictions  could  be  used  to  assess  the  risks  associated  with  increased  development. 

In  order  for  the  goodness-of-fit  approach  to  be  applied,  measurements  and 
predictions  must  be  available  at  the  same  times  and  locations.  If  model  predictions 


FIGURE  3-12 

INPUT  AND  OUTPUT  ERRORS  IN  A  TYPICAL  MODELING  STUDY 


arc  actually  averages  over  extended  spatial  regions,  well  observations  should  be 
averaged  so  that  compatible  regional  measurements  can  be  obtained.  This  may  be 
accomplished  with  a  blocking  scheme,  although  it  is  usually  better  to  smooth  out 
the  measured  head  values  with  a  contour  plot  (sec  Figure  3-11).  The  contoured 
heads  may  then  be  integrated  over  appropriate  model  elements  to  provide  average 
measured  heads. 

Once  a  direct  comparison  of  observations  and  predictions  is  possible,  the 
goodness-of-fit  errors  can  be  computed.  Statistical  error  measures  such  as  the  mean- 
squared-deviation  or  the  value  of  the  Oa^'  percentile  error  may  then  be  derived.  In 
fact,  the  entire  theoretical  framework  of  classical  regression  analysis  may  be  used 
to  analyze  the  “goodness-of-lit”  between  model  and  data. 

The  goodness-of-lit  approach  can  be  useful,  particularly  if  it  is  based  on  a 
careful  regression  analysis,  but  it  has  some  important  practical  limitations.  Most 
of  these  are  related  to  the  fact  that  a  goodness-of-fit  analysis  of  model  performance 
over  a  brief  historical  period  does  not  necessarily  convey  any  information  about  the 
model’s  long-term  prediction  capabilities.  The  extrapolation  of  past  performance 
into  the  future  presumes  that  past  and  future  errors  are  statistically  similar  (i.e., 
arc  drawn  from  the  same  population).  If  the  historical  observation  period  is  brief, 
as  is  usually  the  case  in  groundwater  studies,  it  could  very  well  be  statistically 
anomalous.  If,  for  example,  recharge  and  pumpage  during  the  historical  period 
are  unusually  low,  a  model  with  incorrect  transmissivities  may  give  a  good  fit  to 
measured  heads,  provided  that  its  initial  conditions  and  boundary  conditions  are 
adjusted  appropriately.  The  same  model  inputs  may  give  much  poorer  results  at  a 
later  time  when  pumpage  and  recharge  increase. 

In  order  to  provide  a  reliable  assessment  of  accuracy  the  goodness-of-fit  ap¬ 
proach  requires  a  large  number  of  well  observations.  There  arc  rarely  enough  ob¬ 
servations  available  in  a  typical  groundwater  study  to  support  a  thorough  statistical 
analysis  of  model  fit.  This  is  certainly  the  case  in  the  San  Andres- Glorieta  applica¬ 
tion.  Those  measurements  that  are  available  are  usually  used  to  estimate  aquifer 
parameters  with  a  regional  model  calibration.  If  the  same  measurements  are  also 
used  to  evaluate  accuracy,  they  will  give  a  misleadingly  optimistic  assessment  of  the 
model’s  predictive  ability.  As  an  extreme  example,  it  might  be  noted  that  three  well 
observations  taken  in  different  years  can  be  fit  perfectly  with  a  quadratic  equation 
which  depends  only  on  time.  The  quality  of  this  fit  clearly  docs  not  guarantee 
that  a  quadratic  equation  will  predict  future  water  levels  without  error.  A  much 
better  measure  of  the  model’s  performance  would  be  obtained  if  its  prediction  were 
compared  with  fifty  measurements.  As  a  general  rule,  the  number  of  observations 
used  in  goodness-of-fit  test  should  be  large  compared  to  the  number  of  unknown 
parameters  in  the  model. 

A  promising  alternative  to  the  goodness-of-fit  approach  is  one  based  on  sen¬ 
sitivity  analysis.  Sensitivity  analysis  does  not  require  any  field  data  but  is,  rather, 
based  on  an  investigation  of  the  model’s  response  to  specified  input  errors.  The  basic 
idea  may  be  illustrated  by  considering  a  dynamic  simulation  run  on  a  vertically 
averaged  model  with  one  transmissivity  T  and  one  storage  coefficient  S.  Suppose 
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that  the  predicted  head  at  a  given  time  and  location  is  written  as  h(x,y,t,T,S)  to 
acknowledge  the  model’s  dependence  on  the  inputs  T  and  S.  If  each  of  these  inputs 
were  changed  the  resulting  change  in  the  head  prediction  would  be  approximately: 

Ah(x>y,t,T,S)  =  ^AT+~AS  (3-36) 

Here  the  partial  derivatives  dh/dT  and  dh/dS  are  the  changes  in  head  prediction 
obtained  if  either  the  transmissivity  or  storage  coeflicient  (but  not  both)  is  changed 
infinitesimally.  The  perturbations  AT  and  AS  are  changes  in  transmissivity  and 
storage  coefficient  (not  necessarily  small)  which  shift  the  head  prediction  by  Ah. 

Equation  (3-36)  is  a  linear  (or  first-order)  Taylor  series  approximation  to  the 
model’s  head  solution,  written  as  a  function  of  the  estimated  parameters  T  and  S. 
This  scries  may  be  generalized  to  include  all  model  inputs  as  follows: 

m  dh 

A h(x,y,t)  =  £  ^;Au<  (3-37) 

Here  u,  represents  the  ith  model  input  (out  of  a  total  of  m  inputs)  and  dh/dui  rep¬ 
resents  the  sensitivity  of  the  prediction  at  location  (x,  y)  and  time  t  to  a  small  change 
in  u,.  If  the  input  pertubations  (Atq’s)  are  interpreted  as  errors  (i.c.,  differences 
between  estimated  and  true  values)  then  Ah(x,y,t)  is  the  model’s  prediction  error. 

Equation  (3-37)  may  not  seem  particularly  useful  for  accuracy  analysis  since 
the  input  errors  are  not  known.  This  equation  may,  however,  be  used  to  establish 
an  upper  bound  on  the  prediction  error  if  absolute  values  are  taken  throughout: 

m  3L 

|Ah(z,y,t)|  <  |  — |[e,|  (3-38) 

This  equation  states  that  the  magnitude  of  the  prediction  error  does  not  exceed  the 
summation  on  the  right-hand  side,  with  c,-  selected  as  an  upper  bound  on  the  error 
associated  with  input  u;. 

Equation  (3-38)  may  be  used  to  analyze  model  accuracy  in  the  following  way. 
The  model  is  set  up  and  run  over  the  desired  prediction  period  as  usual.  Then 
its  various  sensitivity  derivatives  are  computed  by  perturbing  each  input  in  turn. 
An  error  bound  is  then  postulated  for  each  model  input.  These  bounds  should 
be  based  on  the  modeler’s  best  estimate  of  the  accuracy  of  the  input  estimate, 
considering  data  availability,  the  degree  of  extrapolation,  and  all  the  other  factors 
mentioned  at  the  beginning  of  this  section.  Equation  (3-38)  is  then  used  to  compute 
the  prediction  error  bound.  The  entire  process  may,  in  principle,  be  automated  and 
a  printout  of  Ah(x,  y,  t)  included  in  the  model’s  output. 

The  sensitivity  approach  has  several  advantages  which  make  it  worth  serious 
consideration.  First,  it  is  not  data-dependent  and  may,  in  fact,  be  used  to  inves¬ 
tigate  accuracy  early  in  the  modeling  study  before  data  are  compiled  and  inputs 
are  estimated.  Second,  it  provides  a  bound  on  prediction  error  which  varies  with 
location  and  time.  This  allows  the  modeler  to  account  for  the  growing  errors  which 
may  occur  when  current  trends  are  projected  far  into  the  future.  Finally,  sensitivity 


analysis  allows  the  modeler  to  identify  the  largest  contributions  to  prediction  error 
(aquifer  parameters,  boundary  conditions,  recharge,  etc.)  so  that  lie  knows  where 
to  focus  his  efforts  to  improve  accuracy. 

The  major  disadvantage  of  sensitivity  analysis  is  its  expense.  The  sensitivity 
derivatives  can  be  quite  expensive  to  compute  if  the  model  has  many  inputs  and  the 
prediction  period  is  long.  The  expense  increases  greatly  if  higher-order  terms  are 
included  in  the  Taylor  series  expansion  since,  in  this  case,  higher-order  sensitivity 
derivatives  must  also  be  derived. 

Limited  sensitivity  analyses  are  beginning  to  appear  in  groundwater  modeling 
studies  (I iconic,  1980;  II  ICC,  1982)  and  there  is  a  growing  recognition  of  the  connec¬ 
tion  between  sensitivity  and  prediction  error.  Hut  systematic  accuracy  evaluations 
of  the  type  outlined  here  arc  still  rare.  At  the  moment,  the  good  ness- of- fit  approach 
remains  the  one  most  accepted. 

3.3.5  Presentation  of  Model  Results 

The  last  phase  of  an  aquifer  water  supply  study  is  the  culmination  of  all  the 
analysis,  interpretation,  and  manipulation  discussed  in  the  preceding  sections.  Here 
the  model’s  predictions  are  generated  and  reported.  If  everything  up  to  this  point 
has  been  done  properly,  the  model  should  yield  a  plausible  set  of  results  which  can 
be  defended  from  available  field  data  and  generally  accepted  theory.  The  modeling 
study  is,  however,  incomplete  if  the  modeler  presents  his  results  in  a  way  which 
implies  that  they  are  perfect.  It  is  equally  incomplete  if  the  modeler  admits  that  his 
results  arc  uncertain  but  does  not  provide  any  further  guidance.  A  good  modeling 
study  acknowledges  and  quantifies  the  uncertainty  of  its  predictions.  This  enables 
potential  critics  to  make  informed  decisions  about  the  model’s  credibility  instead 
of  being  forced  to  either  wholeheartedly  accept  or  completely  reject  its  predictions. 
Such  an  informed  approach  is  ultimately  in  the  best  interest  of  everyone  involved. 


1.  A  COMPARISON  OF  THREE  MODELING  APPROACHES 


The  San  Andres-Glorieta  management  problem  outlined  at  the  beginning  of 
Chapter  2  appears  to  be  a  natural  candidate  for  a  groundwater  modeling  study. 
'Hie  impact  evaluation  question  at  issue  is  clear  but  difficult  to  answer  without  an 
analysis  of  groundwater  How.  Available  aquifer  information  is  not  abundant  but  it 
may  be  adequate  for  a  regional  modeling  study.  Given  this,  it  may  seem  that  any 
experienced  modeler  should  be  able  to  start  with  the  problem  definition  and  data 
base  reviewed  in  Chapter  2,  apply  the  general  principles  of  Chapter  3,  and  end  up 
with  a  set  of  reasonably  credible  predictions.  There  are,  however,  many  judgments 
and  decisions  to  be  made  along  the  way  and  the  predictions  which  emerge  depend 
nearly  as  much  on  the  modeler’s  own  abilities  and  biases  as  on  objective  fact.  This 
is  clearly  illustrated  by  three  modeling  studies  carried  out  to  investigate  the  effects 
of  pumpage  in  the  San  Andres-Cloriela  aquifer.  These  studies  all  relied  on  the 
same  general  principles,  data  base,  and  simulation  model  but  yielded  very  different 
predictions;  so  different  that  they  virtually  cancelled  each  other  out. 

The  simulation  approaches  used  in  the  three  San  Andres-Glorieta  studies  are 
distinguished  primarily  by  differences  in  model  formulation  and  input  estimation. 
A  brief  summary  of  each  approach  is  provided  below: 

Modeling  Study  1  (Geohydrology, 1982) 

This  study  solved  the  drawdown  equation  directly.  No-flow  boundaries  were 
assumed  on  all  sides  of  the  solution  region.  The  aquifer  parameters  were 
either  postulated  or  estimated  from  pump  test  data.  Incremental  pumpage 
was  obtained  from  the  Plains  Electric  pumping  schedule  (see  Figure  2-2)  and 
incremental  recharge  was  assumed  to  be  zero. 

Modeling  Study  2  (Geotrans,1982) 

Drawdown  was  computed  in  this  study  from  a  simulation  initialized  with 
steady-state  heads  based  on  1968  pumpage  and  recharge  values.  Using  the 
terminology  of  Section  3.2.5,  the  nominal  pumping  strategy  was  defined  by 
1968  conditions  and  the  modified  strategy  by  1968  conditions  with  Plains 
Electric  pumpage  added.  No-flow  boundaries  were  assumed  on  all  sides  of 
the  solution  region.  Aquifer  parameters  were  estimated  from  pump  test  data 
and  a  steady-state  regional  calibration. 

Modeling  Study  3  (HEC,1982) 

This  study  used  the  same  differencing  approach  as  Modeling  Study  2  to  derive 
drawdown  predictions.  Specified  head  boundaries  were  used  on  three  sides 
of  the  solution  region  and  a  no- flow  boundary  on  the  fourth  side.  Aquifer 
parameters  and  boundary  heads  were  either  postulated  or  derived  from  a 
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steady-state  calibration.  This  study  differed  from  the  other  two  in  that  the 
emphasis  was  on  the  impact  of  pumpage  on  water  levels  at  Fort  Wingate  Army 
Depot  (located  approximately  10  miles  west  of  the  Ciniza  well  field).  This  had 
some  influence  on  the  way  the  model  was  formulated  and  on  the  selection  of 
model  input  values. 

The  maximum  drawdowns  predicted  by  each  model  are  plotted  in  Figure  4-1. 
Note  that  the  area  pictured  in  this  figure  covers  the  south-central  portion  of  the 
San  Andres-Gloricta  region.  Each  of  the  three  models  predicted  that  maximum 
drawdown  will  occur  after  about  18  years  of  plant  operation.  It  is  apparent  that 
Model  2  predicts  the  greatest  impact  (the  1000  foot  drawdown  indicated  at  the  well 
field  locally  dewaters  the  aquifer).  Model  3  predicts  relatively  modest  and  localized 
drawdowns  while  Model  1  gives  results  somewhere  in  between.  The  maximum 
difference  between  these  predictions  in  the  vicinity  of  the  well  field  is  over  600  feet. 

The  analysis  presented  in  this  chapter  attempts  to  discover  how  three  modeling 
studies  using  such  similar  methodologies  can  generate  such  different  results.  The 
case  study  comparison  is  interesting  for  its  own  sake  but  it  also  reveals  something 
about  the  practice  of  groundwater  modeling  in  general.  The  data  ambiguities  and 
opportunities  for  subjective  interpretation  encountered  in  the  San  Andres-Glorieta 
play  an  equally  important  role  in  other  water  supply  investigations.  The  differences 
examined  here  are  dramatic  but  not  unusual. 

The  comparison  begins  with  a  description  of  the  general  specifications  for 
each  model — the  simulation  approach  used,  the  boundaries  of  the  solution  region, 
etc.  It  then  considers  the  important  application  issues  addressed  in  Section  3.3— 
discretization,  input  estimation,  and  accuracy  evaluation.  The  chapter  concludes 
with  a  review  of  model  similarities  and  differences  and  an  overall  assessment  of  the 
case  study. 

4.1  MODEL  FORMULATION 

All  three  of  the  studies  reviewed  here  adopted  the  two-dimensional  vertically 
averaged  modeling  approach  described  in  Section  3.2.  In  each  case,  the  flow  equa¬ 
tion  was  solved  with  the  USGS  two-dimensional  finite  difference  program  docu¬ 
mented  in  Trcscott,  Pinder,  and  Larson  (1976).  This  program  can  handle  both 
confined  and  unconfined  conditions  and  provides  reasonable  flexibility  for  laying 
out  boundaries  and  assigning  boundary  conditions. 

The  solution  regions  for  the  three  models  represent  somewhat  different  ap¬ 
proaches  to  aquifer  simulation  (see  Figure  4-2).  The  network  for  Model  1  is  essen¬ 
tially  rectangular  with  some  modification  along  the  lower  (southwestern)  boundary 
to  account  for  the  irregular  nature  of  the  aquifer  outcrop.  The  areal  extent  of  this 
network  covers  much  of  the  aquifer  but  does  not  actually  reach  generally  acknowl¬ 
edged  geological  boundaries  except  along  the  southern  edge.  The  limited  coverage 
of  the  model  network  can  be  justified  if  the  fluxes  crossing  the  right  (eastern),  left 


C)  model  number  Z> 
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MAXIMUM  PREDICTED  DRAWDOWNS  FOR  MODELS  1,2.  AND  3 


(western),  or  upper  (northern)  boundaries  arc  reasonably  constant  and  arc  unaffected 
by  purnpage  at  the  Ciniza  well  field.  Model  1  assumes  that  this  is  the  ease  and, 
consequently,  imposes  zero  gradient  drawdown  conditions  on  all  four  sides  of  the 
network. 

The  nearly  rectangular  network  used  in  Modeling  Study  2  extends  north  to 
the  edge  of  the  aquifer  but  approximates  the  lower  boundary  somewhat  differently 
than  Model  1.  Model  2  assumes  that  there  is  no  flow  across  any  of  the  network 
boundaries. 

The  perfectly  rectangular  Model  3  network  makes  no  attempt  to  cover  the 
entire  aquifer  but  includes  only  the  region  where  drawdowns  are  expected  to  be 
significant.  The  influence  of  the  surrounding  flow  field  is  accounted  for  by  specified 
head  boundary  conditions  on  the  upper,  lower,  and  left  sides  of  the  network.  The 
right  side  is  assumed  to  be  a  no- flow  boundary.  The  boundary  heads  of  Model  3 
are  adjusted  so  that  groundwater  enters  the  network  from  the  outcrop  area,  moves 
north,  and  gradually  curves  to  the  west.  This  assumed  flow  pattern  is  based  on 
water  level  data  reported  by  Shomakcr  (1971)  and  plotted  in  Figure  2-5.  Since  the 
boundary  heads  are  assigned  constant  values,  they  are  implicitly  assumed  to  be 
unaffected  by  purnpage  at  Ciniza. 

The  zero-gradient  boundary  condition  for  Model  1  is  a  straightforward  con¬ 
sequence  of  the  drawdown  formulation  outlined  in  Section  3.2.5.  The  boundary 
conditions  for  Models  2  and  3  arc  somewhat  more  subjective  and  have  broader  im¬ 
plications.  The  zero-flux  approach  taken  in  Model  2  treats  the  aquifer  as  a  closed 
system  which  interacts  with  the  outside  world  only  through  the  horizontal  boundary 
flux  (called  qu  in  Section  3.2.3).  In  this  case  the  aquifer  can  be  in  steady-state  only 
if  the  purnpage,  recharge,  and  leakage  components  of  the  horizontal  flux  balance 
exactly.  Since  purnpage  and  recharge  are  specified  independently,  the  leakage  term 
must  adjrst  to  give  a  steady-state  solution.  The  open  (specified  head)  boundaries 
used  in  Model  3  provide  a  different  way  to  achieve  a  steady-state  balance.  In  this 
case  the  head  gradients  along  each  of  the  specified  head  boundaries  may  adjust  to 
allow  groundwater  to  either  enter  or  leave  the  network.  The  difference  between 
purnpage  and  recharge  does  not  have  to  be  absorbed  by  leakage  across  the  horizon¬ 
tal  boundary.  The  leakage  term  of  Model  2  and  the  open  boundaries  of  Model  3 
are  both  reasonable  ways  to  account  for  interactions  between  the  aquifer  and  its 
environment.  They  are,  however,  clearly  not  equivalent. 

A  comparison  of  the  three  model  formulations  suggests  that  Model  2  takes  the 
safest  (most  conservative)  approach  since  it  makes  the  fewest  assumptions.  Both 
Model  1  and  Model  3  require  the  network  boundaries  to  be  unaffected  by  purnpage 
at  the  Ciniza  well  field.  Model  3  also  requires  estimation  of  head  boundary  values 
along  three  sides  of  the  solution  region.  Although  such  requirements  may  not 
necessarily  pose  problems,  it  is  best  to  avoid  them  if  possible. 
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4.2  MODEL  APPLICATION 


4.2.1  Spatial  and  Temporal  Discretization 

Since  the  three  case  study  models  were  all  run  with  the  same  finite  difference 
computer  program,  they  all  used  the  same  discretization  approach.  Each  solution 
region  was  divided  into  a  grid  of  parallel  lines  spaced  close  together  near  the  middle 
of  the  region  and  progressively  further  apart  toward  the  edges  (see  Figure  3-8b  for 
an  illustration  of  this  approach).  The  smallest  and  largest  grid  intervals  for  each 
model  are  listed  below: 


MODEL 

SMALLEST  INTERVAL 

LARGEST  INTERVAL 

(miles) 

(miles) 

1 

0.25 

10.0 

2 

1.00 

10.0 

3 

0.25 

2.0 

The  smallest  interval  defines  the  model’s  spatial  resolution  near  the  Ciniza  well 
field.  Temporal  discretization  for  each  model  followed  the  approach  illustrated  in 
Figure  3-9.  Model  inputs  are  assumed  constant  over  each  time  step  and  model 
predictions  are  assumed  to  vary  linearly.  The  model  time  steps  were  all  one  year 
or  less.  Since  the  time  horizon  of  interest  in  the  impact  evaluation  is  forty  years, 
annual  time  steps  probably  provide  adequate  temporal  resolution. 

The  discretization  approach  used  in  each  of  the  three  models  is  reasonable  and 
adequate  for  impact  evaluation  purposes.  It  seems  safe  to  say  that  discretization  had 
no  significant  effect  on  the  differences  in  model  predictions  which  are  of  particular 
interest  here. 

4.2.2  Input  Estimation 

Given  the  general  similarity  of  the  model  formulations  and  discretization  pro¬ 
cedures  outlined  above,  it  seems  obvious  that  the  divergent  predictions  generated  by 
the  three  case  study  models  were  caused  by  differences  in  model  inputs,  i.e.  in  the 
values  used  for  aquifer  parameters,  boundary  conditions,  and  pumpage/recharge 
rates.  What  is  less  obvious  is  why  qualified  hydrologists  should  derive  such  different 
input  estimates  from  the  same  data  base.  This  issue  is  clearly  worth  further  inves¬ 
tigation. 

It  is  helpful  to  summarize  the  inputs  that  need  to  be  estimated  for  each 
model  since  differences  in  formulation  lead  to  differing  input  requirements.  These 
requirements  are  indicated  by  closed  circles  in  the  following  table: 


en 


INPUT 


1 


MODEL 

2 


3 


Aquifer  parameters  • 

Initial  heads 
Adjoining  aquifer  head 
Boundary  heads 


Nominal  pumpage  •  • 

Nominal  recharge  • 


In  order  to  keep  the  differences  between  the  three  modeling  studies  clear,  it  is 
convenient  to  examine  their  input  estimation  methods  separately.  Each  model  is 
covered  in  one  of  the  subsections  which  follow. 


Model  1 

Aquifer  parameter  estimates  for  Model  1  were  obtained  from  a  variety  of 
primary  and  secondary  sources  which  are  not  very  well  documented.  Transmissivities 
were  assumed  to  be  uniform  over  large  blocks  as  shown  in  Figure  4-3.  The  trans¬ 
missivity  estimates  used  vary  dramatically  from  a  low  of  70  ft2 /day  to  over  77,000 
ft2/day.  These  estimates  appear  to  have  been  extrapolated  from  individual  pump 
test  or  soil  sample  values  by  means  of  a  blocking  scheme  similar  to  the  one  il¬ 
lustrated  in  Figure  3- lib. 

The  procedure  used  to  estimate  the  transmissivity  value  for  the  Ciniza  area  is 
worth  considering  in  detail  since  this  value  has  a  significant  effect  on  the  model’s 
drawdown  predictions.  Since  1956  Shell  Oil  has  pumped  about  600  acre-feet/year 
from  a  well  field  near  Ciniza.  Water  level  observations  collected  at  several  nearby 
wells  over  the  1956-1982  period  indicate  a  gradual  decline  in  head  which  is  presum¬ 
ably  due  to  the  Shell  pumpage.  If  the  1956  water  level  is  taken  as  a  reference, 
drawdowns  for  the  historical  period  may  be  computed  and  plotted  vs.  time.  A 
conventional  pump  test  analysis  of  this  plot  carried  out  in  Modeling  Study  1  gave 
a  transmissivity  estimate  of  about  450  ft2/day.  This  was  apparently  confirmed  by 
a  limited  regional  calibration  study. 


FIGURE  4-3 

TRANSMISSIVITY  REGIONS  FOR  MODEL  1 


It  should  be  noted  that  other  pump  test  analyses  of  the  Ciniza  drawdown  data 
produced  transmissivity  estimates  which  varied  from  less  than  80  ft“/day  to  over 
3000  ft“/day  (Dames  and  Moore,  1982).  The  wide  range  of  uncertainty  observed 
in  these  estimates  reflects  different  assumptions  made  in  type  curve  computations. 
The  average  difference  in  computed  transmissivity  (about  1500  ft2/day)  could  be 
taken  as  an  upper  bound  on  the  transmissivity  estimation  error  for  the  Ciniza  area. 

A  confined  storage  coefficient  value  of  5  X  10-'1  (unitlcss)  for  the  Ciniza  area 
was  derived  from  the  Shell  Oil  pump  test  analysis  mentioned  earlier.  This  value 
was  used  throughout  the  confined  portion  of  the  aquifer  but  was  increased  to  0.05 
in  the  unconCncd  outcrop  area.1  This  increase  accounts  for  the  impact  which  a 
rising  or  falling  water  table  has  on  storage.  It  has  the  effect  of  making  local  heads 
less  responsive  to  changes  in  pumpage.  The  0.05  estimate  for  the  outcrop  storage 
coefficient  is  plausible,  considering  that  specific  yields  for  unconfincd  sandstone  and 
limestone  formations  usually  fall  in  the  range  from  0.02  to  0.08  (Linsley,  et  al., 1982). 

The  Model  1  lcakance  value  was  computed  by  assuming  that  the  confining  layer 
lying  between  the  San  Andres-Gloricta  and  the  adjoining  aquifer  (presumably  the 
Sonsela  formation)  has  a  vertical  hydraulic  conductivity  of  1  X  10~l2ft/sec.  and 
is  600  feet  thick.  This  gives  a  leakance  estimate  of  1.7  X  10~15  sec-1  or  about 
1.4  X  10~10  days-1  (see  Equation  3-11).  The  computed  leakance  was  assumed  to 
apply  uniformly  throughout  the  Model  1  network. 

It  is  instructive  to  see  what  this  leakance  value  implies  about  the  role  of  leakage 
in  the  Model  1  water  budget.  The  incremental  decrease  in  leakage  out  of  the  aquifer 
due  to  the  Plains  Electric  pumpage  may  be  calculated  as  follows: 

A q,  =  LAd  (4  -  1) 

where  L  is  the  lcakance  and  d  is  the  average  drawdown  over  an  area  A.  The 
maximum  Model  1  drawdown  prediction  plotted  in  Figure  4-1  can  be  roughly 
approximated  as  250  feet  over  a  circle  with  a  radius  of  12  miles.  Given  this 
approximation  and  a  leakance  value  of  1.4  X  10~10  days-1,  Equation  (4-1)  yields 
a  decrease  in  leakage  of  3.7  acre-feet/year.  When  compared  to  an  incremental 
pumpage  value  of  4000  acre-feet/year  this  is  negligible.  Leakage  clearly  has  little 
effect  on  the  predictions  of  Model  1. 


*Thc  Model  1  outcrop  area  covers  the  two  transmissivity  blocks  at  the  bottom  of  Figure  4-3  labeled 
7754  and  43,817. 
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Model  2 


Aquifer  parameters  for  Model  2  were  estimated  primarily  from  pump  test 
results  and  qualitative  geological  information.  In  some  eases,  these  estimates  were 
refined  with  a  regional  model  calibration  based  on  steady-state  head  observations. 
The  Model  2  aquifer  parameters  were  extrapolated  with  a  contouring  procedure 
which  allows  the  parameters  to  vary  continuously  throughout  the  aquifer,  as  is 
illustrated  by  the  transmissivity  map  of  Figure  4-4. 

A  comparison  of  Figures  4-3  and  4-4  reveals  how  different  modelers  can  inter¬ 
pret  and  extrapolate  the  same  set  of  pump  test  results.  The  two  transmissivity  maps 
arc  qualitatively  similar — in  each  case  transmissivities  around  the  top  and  two  sides 
of  the  aquifer  are  low  and  transmissivities  near  the  eastern  end  of  the  outcrop  are 
very  high.  The  most  important  interpretive  difference  is  in  the  Ciniza  area,  where 
the  Model  2  contours  bend  in  such  a  way  as  to  lower  the  transmissivity  north  of 
the  well  field  from  the  Model  1  value  of  450  ft2/day  to  less  than  100  ft2/day.  The 
low  transmissivities  used  in  Model  2  were  obtained  from  an  independent  analysis 
of  pump  test  data  from  the  Ciniza  well  field. 

Confined  storage  coefficients  for  Model  2  were  computed  by  multiplying  a 
specific  storage  of  4  X  10~7  ft-1  by  the  estimated  aquifer  thickness  (in  ft.).  This 
gives  values  which  vary  throughout  the  aquifer  but  are  generally  less  than  1  X  10~4 
near  the  Ciniza  well  field.  The  4  X  10-7ft_l  value  used  for  specific  storage  seems 
rather  low.  A  rough  check  on  this  value  can  be  obtained  by  computing  the  specific 
storage  from  Equation  (2-5).  The  range  for  the  bulk  compressibility  of  a  jointed 
rock  aquifer  given  by  Freeze  and  Cherry  (1979)  is  10”8  to  10~10  m2/Newton  (about 
lfT12  to  10~u  psi_l).  If  the  smaller  value  is  used  with  a  porosity  of  0.10,  the 
resulting  specific  storage  is  1.5  X  10_®  ft-1.  This  compares  well  with  the  value  of 
1.0  X  10-6ft-1  given  by  Lohman  (1972).  If  1.5  X  10_6ft_1  is  used,  the  computed 
storage  coefficient  for  the  Ciniza  area  quadruples  to  4  X  10-4.  This  is  close  to 
the  value  used  in  Model  1.  The  unconfined  storage  coefficient  for  Model  2  was 
assumed  to  have  a  uniform  value  of  0.10  wherever  unconfincd  conditions  exist.  This 
is  consistent  with  the  range  of  reasonable  values  cited  in  the  Model  1  discussion. 

The  leakance  value  used  for  most  of  the  Model  2  network  was  rather  arbitrarily 
set  equal  to  1.0  X  10~,4sec_1.  Sensitivity  runs  showed  that  this  small  leakance 
value  has  little  effect  on  the  model’s  predictions.  The  leakance  must  be  increased 
to  roughly  1.0  X  10-,3sec-1  before  the  resulting  decrease  in  leakage  flux  becomes 
significant  compared  to  the  Plains  Electric  pumping  rate.  Model  2  used  a  much 
higher  leakance  value  (1.0  X  10~,0sec"1)  along  the  Nutria  monocline  (a  structural 
feature  located  about  four  miles  east  of  Gallup)  but  the  resulting  effect  on  the 
predicted  head  is  very  localized.  For  all  practical  purposes,  leakage  plays  a  negligible 
role  in  the  Model  2  analysis. 

Since  Model  2  assumes  no-flow  boundaries  on  all  four  sides  the  only  inputs 
required  for  the  steady-state  head  simulation  (other  than  aquifer  parameters)  are 
the  nominal  recharge  and  pumpage  and  the  head  in  the  adjoining  aquifer  (the 
Sonsela  formation).  As  was  pointed  out  in  Section  4.1,  an  aquifer  surrounded  by 


no-flow  boundaries  can  be  in  steady-state  only  if  leakage  is  able  to  account  for  the 
difference  between  punipage  and  recharge.  A  simple  water  balance  (sec  Table  2-2) 
reveals  that  San  Andre9-Gloricta  recharge  exceeded  pumpage  in  1968  by  several 
thousand  acre-feet.  Since  the  Model  2  leakancc  values  are  small,  the  head  in  the 
Sonsela  would  have  to  be  at  least  several  hundred  feet  lower  than  the  head  in  the 
San  Andres-Gloricta  in  order  for  this  much  water  to  leave  by  leakage.  This  is 
unrealistic  considering  that  observed  head  differences  between  the  Sonsela  and  San 
Andres-Gloricta  are  only  about  150  feet.  This  dilemma  was  resolved  by  inserting 
an  ad  hoc  “line  of  discharge”  north  of  the  outcrop  area.  SpcciGed  discharge  values 
along  this  line  force  enough  water  out  of  the  San  Andres-Gloricta  to  achieve  a 
steady-state  consistent  with  the  150  foot  observed  head  difference  mentioned  above. 
Additional  lines  of  discharge  were  also  apparently  placed  below  eastern  sections  of 
Rio  San  Jose. 

The  nominal  recharge  inputs  used  in  Model  2  were  based  on  Shomaker’s  (1971) 
estimates  of  average  infiltration  in  the  outcrop  area.  A  fixed  infiltration  rate  of  0.75 
inches/ycar  was  assumed  to  apply  over  an  outcrop  area  covering  approximately 
250  square  miles.  This  gives  a  total  recharge  flux  of  about  10,000  acre-feet/year, 
as  compared  to  the  value  of  8040  acre-feet/year  reported  in  Table  2-2.  Nominal 
pumpage  rates  assumed  in  the  Model  2  steady-state  simulation  were  obtained  from 
generally  available  records  and  are  consistent  with  the  values  presented  in  Table 
2-1. 

The  steady-state  simulation  results  obtained  with  the  discharge  lines,  recharge 
rates,  and  purnpages  assumed  in  Model  2  are  shown  in  Figure  4-5.  The  Model  2 
simulation  implies  that  groundwater  moves  from  the  western  end  of  the  outcrop 
northward  (as  suggested  by  Shomaker,  1971)  and  then  swings  southward  toward 
Rio  San  Jose.  This  flow  pattern  is  credible  in  the  area  where  head  observations  are 
available  but  it  is  difficult  to  accept  in  the  east  where  the  flow  reverses  direction. 
It  appears  that  the  convenient  but  speculative  lines  of  discharge  dominate  the 
regional  flow  field,  forcing  the  suprising  reversal.  Although  the  steady-state  head 
distribution  simulated  by  Model  2  may  be  correct,  it  depends  on  some  rather 
arbitrary  assumptions.  The  “line  of  discharge”  approach  needs  to  be  justified 
by  geological  and  hydrological  data  since  it  could  clearly  be  abused  if  applied 
indiscriminately.  In  this  case  the  “line  of  discharge”  seems  to  be  primarily  a  method 
for  resolving  the  problems  which  result  from  the  imposition  of  universal  no-flux 
boundary  conditions  on  a  steady-state  model.  This  is  clearly  not  an  adequate 
justification,  since  these  problems  can  be  best  handled  by  adopting  more  realistic 
(and  more  flexible)  boundary  conditions. 
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FIGURE  4-5 

STEADY-STATE  SIMULATION  RESULTS  FOR  MODEL  2 


Model  3 

Aquifer  parameters  for  Model  3  were  initially  obtained  from  an  early  version  of 
Model  1  and  then  were  refined  with  a  regional  model  calibration  based  on  the  1968 
head  observations  discussed  earlier.  The  Model  3  transmissivities  were  blocked  in 
a  manner  similar  to  those  of  Model  1  (see  Figure  4-6).  The  most  notable  feature 
of  the  Model  3  transmissivity  distribution  is  the  higher  value  used  in  the  Ciniza 
vicinity  (1500  ft2/day).  This  is  at  the  upper  end  of  the  range  of  estimates  derived 
from  the  Ciniza  pump  test  data. 

Model  3  assumes  a  value  of  5  X  10-4  for  the  confined  storage  coefficient  (the 
same  as  Model  1).  An  unconfined  value  is  not  required  since  the  Model  3  network 
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docs  not  extend  into  tlic  outcrop  area.  It  should  be  noted  that  Model  3  uses 
an  open  (specified  head)  southwestern  boundary  to  account  for  inflows  from  the 
outcrop  region.  These  inflows  automatically  adjust  to  balance  outflows  across  other 
boundaries  and  withdrawals  due  to  pumpage.  Since  this  approach  is  intended  to 
account  for  both  leakage  and  recharge,  the  leakance  coefficient  and  nominal  recharge 
rate  for  Model  3  were  assumed  to  be  zero.  Overall,  the  Model  3  aquifer  parameters 
arc  as  reasonable  as  those  of  the  other  models,  with  the  notable  exception  of  the 
Ciniza  transmissivity,  which  may  be  as  much  as  an  order  of  magnitude  too  high. 
The  implications  of  this  assumption  are  considered  in  Section  4.3. 

4.2.3  Accuracy  Evaluation 

The  approach  to  accuracy  evaluation  taken  in  Modeling  Studies  1  and  2  rep¬ 
resents  what  might  be  called  the  “model  validation”  point-of-view.  This  approach 
is  concerned  primarily  with  proving  that  the  model  is  an  acceptable  description 
of  reality.  Once  the  proof  of  validity  is  accepted  the  predictions  are  presumed  to 
be  sufficiently  accurate  for  their  intended  purpose.  This  all-or-nothing  approach 
provides  no  way  to  objectively  compare  different  models.  All  validated  models  are, 
by  implication,  equally  accurate. 

The  goodness-of-fit  and  sensitivity  analysis  methods  described  in  Section  3.3.4 
represent  a  more  realistic  and  informative  approach  to  the  accuracy  issue.  The  goal 
of  these  methods  is  to  derive  upper  bounds  on  the  model’s  prediction  errors.  When 
expressed  as  functions  of  time  and  space  these  bounds  clearly  reveal  the  model’s 
strengths  and  weaknesses. 

Since  there  are  not  enough  head  measurements  in  the  San  Andres-Glorieta 
aquifer  to  support  an  adequate  goodness-of-fit  investigation,  sensitivity  analysis 
must  be  used  to  evaluate  model  accuracy.  The  “decision  tree”  search  conducted 
in  Modeling  Study  3  is  a  type  of  sensitivity  analysis  which  is  particularly  useful 
in  impact  assessments.  The  decision  tree  approach  attempts  to  identify  the  set 
of  plausible  inputs  which  gives  the  worst  case  impact.  In  Modeling  Study  3  trans¬ 
missivities,  storage  coefficients,  and  boundary  conditions  were  progressively  varied 
above  and  below  their  original  values  and  at  each  stage  the  value  giving  the  greatest 
drawdown  was  retained.  The  difference  between  the  model’s  original  drawdown 
prediction  and  the  worst  case  result  (evaluated  at  Fort  Wingate)  was  about  160  feet. 
This  can  be  viewed  as  an  upper  bound  on  the  prediction  error.  It  should  be  pointed 
out  that  the  choice  of  nominal,  upper,  and  lower  values  for  each  parameter  in  a 
sensitivity  analysis  is  usually  rather  subjective.  One  investigator’s  worst  case  may 
be  another’s  nominal  or  best  case.  The  worst  case  values  of  the  Model  3  sensitivity 
analysis  (a  transmissivity  of  250  ft2/day  and  a  storage  coefficient  of  5  X  10~5)  are, 
for  example,  close  to  the  nominal  values  used  in  Model  2.  The  Model  2  investigators 
presumably  feel  these  values  are  reasonable  (“most  likely”)  estimates  rather  than 
worst  case  extremes.  Such  differences  in  interpretation  do  not  discredit  sensitivity 
analysis  as  a  modeling  technique  but  they  do  suggest  that  sensitivity  results  (like 
all  modeling  results)  should  be  carefully  examined  before  they  are  used  to  guide 
policy  decisions. 


4.3  SUMMARY  AND  OVERALL  ASSESSMENT 


When  all  of  the  assumptions,  interpretations,  and  results  of  the  San  Andres- 
Gloricta  case  study  are  reviewed,  it  becomes  apparent  that  many  of  the  most 
noticeable  differences  between  the  three  models  are  not  very  important.  The 
shape,  resolution,  and  extent  of  the  model  network  have  little  impact  since  the 
drawdown  cone  is  contained  well  within  the  boundaries  of  all  three  alternatives. 
The  analysis  presented  in  Section  3.2.4  indicates  that  only  the  type  of  boundary 
condition  (drawdown  vs.  drawdown  gradient)  influences  drawdown  predictions. 
Specific  boundary  values  do  have  an  effect  on  the  simulated  head  which  may,  in  turn, 
influence  the  transmissivity  estimates  generated  in  a  regional  model  calibration.  But 
regional  calibration  was  not  relied  upon  extensively  in  the  modeling  studies  reviewed 
here.  It  was  not  used  at  all  in  Study  1  and  had  relatively  little  impact  on  estimated 
well  field  transmissivities  in  Studies  2  and  3.  All  three  studies  were  designed, 
intentionally,  to  minimize  the  impact  of  boundary  conditions  on  the  drawdown 
results  of  most  interest  from  a  policy  point-of-view.  Similar  comments  apply  to 
the  nominal  pumpage  and  recharge  values  used  to  simulate  existing  hydrologic 
conditions.  Here  again,  these  inputs  have  an  effect  on  head  but  not  on  drawdown. 


It  follows  from  the  brief  review  presented  above  that  the  dominant  inputs  in 
the  San  Andres-Glorieta  case  study  must  be  the  aquifer  parameters,  particularly  the 
transmissivity  and  storage  coefficient.  This  conclusion  should  come  as  no  suprise 
since  it  was  clearly  anticipated  in  Section  3.2.4.  It  is  possible,  however,  to  go 
still  further  and  state  that  only  the  transmissivity  and  storage  coefficient  values 
in  the  vicinity  of  the  Ciniza  well  field  have  a  significant  impact  on  drawdown. 
This  assertion  is  reinforced  by  the  simple  “desk  top”  analysis  summarized  in  the 
following  paragraphs. 

The  importance  of  the  aquifer  parameters  in  the  Ciniza  area  is  related  to  the 
nature  of  the  pumping  strategy  investigated  in  the  case  study.  The  pumping  wells 
to  be  used  by  Plains  Electric  are  clustered  in  a  small  region  located  well  inside 
the  boundaries  of  the  aquifer  in  a  relatively  geologically  homogeneous  region.  It  is 
interesting  to  note  that  the  predicted  drawdown  contours  plotted  in  Figure  4-1  are 
nearly  concentric  circles  centered  on  the  Plains  well  field.  This  suggests  that  the 
impact  evaluation  problem  can  be  solved  with  the  classical  well  drawdown  radial 
flow  equation  (Freeze  and  Cherry,  1979).  Various  solutions  to  this  equation  are 
available  but  it  is  most  instructive  to  begin  with  the  simple  approximation  proposed 
by  Jacob(l940): 


,  2.3A Q ,  2.257T 

d  —  log 


(4-2) 


4* T  r2S 

Here  T  and  S  are  the  transmissivity  and  storage  coefficient  near  a  well  pumping  at 
a  volumetric  rate  AQ.  The  drawdown  d  is  evaluated  at  time  t  at  a  distance  r  from 
the  well.  All  variables  are  assumed  to  be  in  consistent  units. 


Jacob’s  solution  clearly  reveals  how  drawdown  depends  on  the  aquifer  para¬ 
meters  T  and  S.  It  is  also  consistent  with  the  results  of  the  case  study — the 
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TABLE  4-1 


COMPARISON  OF  SIMULATED  DRAWDOWNS  WITH  VALUES  OBTAINED 

FROM  A  THE  IS  SOLUTION 

All  drawdowns  computed  at  t  =  18  years 
Incremental  puinpage  assumed  to  be  5000  acre-fect/year 


predicted  drawdown  was  greater  when  T  and  S  were  assumed  to  be  small  (Model 
2)  and  less  when  T  and  S  were  assumed  to  be  large  (Model  3).  This  is  a  physically 
plausible  result  since  head  gradients  must  become  steeper  as  conductivity  and 
compressibility  decrease  if  a  given  pumping  rate  is  to  be  sustained. 

Jacob’s  approximation  requires  the  argument  of  the  log  function  to  be  small 
and  so  is  not  accurate  if  Tt  is  large  compared  to  r2S.  A  more  general  but  less 
descriptive  well  drawdown  solution  has  been  proposed  by  Thcis  (1935).  The  Theis 
solution  may  be  compared  directly  with  the  model  predictions  plotted  in  Figure  4-1 
if  t  is  set  equal  to  the  time  of  maximum  drawdown  (18  years)  and  r  is  varied  over  a 
range  of  appropriate  distances  (e.g.  4,  8,  and  12  miles).  Each  of  the  three  modeling 
alternatives  may  be  characterized  by  the  combination  of  spatially  averaged  T  and  S 
values  listed  in  Table  4-1.  The  simulated  and  Theis  drawdown  predictions  displayed 
in  this  table  are  suprisingly  close,  confirming  the  hypothesis  that  the  Plains  Electric 
problem  is  a  simple  example  of  drawdown  from  a  pumping  well. 

It  might  be  argued  that  the  only  way  to  know  that  a  radial  flow  solution  is 
appropriate  for  this  problem  is  to  perform  a  complete  two-dimensional  simulation 
first.  That  may  be,  but  it  seems  likely  that  a  Theis  solution  would  occur  to  many 
groundwater  hydrologists  presented  with  the  problem  description  summarized  in 
Chapter  2.  It  is  probably  more  accurate  to  say  that  a  desk- top  analysis  based  on 
the  Theis  solution  simply  isn’t  as  impressive  as  a  full-blown  computer  simulation. 
It  is  sometimes  difficult  to  resist  the  belief  that  a  computer  model  is  better  just 
because  it  is  more  complicated. 

Certainly  there  arc  many  situations  where  computer  models  are  the  only 
reasonable  way  to  analyze  a  complex  problem.  One  case  study  does  not  prove 
that  modeling  is  unnecessary  or  redundant.  Nevertheless,  the  San  Andrcs-Glorieta 
experience  illustrates  the  value  of  carefully  thinking  about  a  problem  before  un¬ 
dertaking  an  expensive  modeling  project.  A  few  simple  hand  calculations  based 
on  readily  available  data  can  help  to  focus  attention  on  the  truly  crucial  aspects 
of  a  problem  and  perhaps  save  months  of  work.  In  the  San  Andres- Glorieta  case 
study,  the  dramatic  differences  in  predicted  drawdowns  were  ultimately  the  result  of 
different  interpretations  of  a  handful  of  pump  tests.  The  accuracy  of  these  predic¬ 
tions  would  probably  have  improved  significantly  if  more  effort  had  been  spent  on 
field  studies  in  the  Ciniza  area. 

In  the  last  analysis,  the  critical  issue  in  this  case  study  was  not  model  for¬ 
mulation  but  data  collection  and  interpretation.  Modeling  has  made  groundwater 
hydrology  more  technical  and  sophisticated  but  it  has  not  eliminated  the  need  for 
carefully  designed  and  executed  field  studies.  Ideally,  the  modeler  and  field  inves¬ 
tigator  should  work  together  so  that  the  strengths  of  their  two  disciplines  can  be 
combined.  Groundwater  models  can  be  used  to  help  design  field  sampling  programs 
by  showing  where  additional  data  will  be  most  beneficial.  Conversely,  sampling 
programs  tailored  to  the  needs  and  assumptions  of  discretized  regional  models  can 
provide  more  reliable  information  for  estimating  model  inputs.  Such  a  cooperative 
approach  is  ultimately  the  most  reasonable  and  cost-effective  way  to  obtain  credible 
impact  predictions. 


5.  GUIDELINES  FOR  GROUNDWATER  MODELING 


The  preceding  chapters  include  a  number  of  recommendations  and  guideline 
which  relate  to  the  practical  details  of  model  application.  These  cover  topics  ranging 
from  the  specification  of  boundary  conditions  to  the  selection  of  an  extrapolation 
technique.  This  chapter  attempts  to  step  back  from  technical  details  and  look  at 
some  of  the  major  decisions  to  be  made  in  a  groundwater  modeling  study.  The 
guidelines  which  follow  are  inspired  by  the  San  Andres-Glorieta  case  study  but 
are  meant  to  be  generally  applicable.  Each  modeler  must,  of  course,  adapt  these 
general  guidelines  to  the  unique  requirements  of  his  particular  problem.  This  could, 
in  fact,  be  viewed  as  the  “first  guideline”. 

The  major  recommendations  of  this  report  can  be  summarized  as  follows: 

1.  Clearly  define  the  objectives  of  the  modeling  study  at  the  outset.  If  the  study 
results  are  to  be  used  for  hydrologic  impact  evaluation  determine  what  type  of 
predictions  are  needed  (head,  drawdown,  storage  changes,  etc.) 

2.  Consult  closely  with  a  groundwater  geologist  familiar  with  the  region  of  inter¬ 
est.  Local  geological  experience  can  help  make  data  interpretation  and  input 
estimation  more  informed  and  realistic. 

3  Always  carry  out  a  thorough  desk- top  analysis  before  starting  a  modeling 
project.  In  particular: 

•  Construct  a  rough  aquifer-scale  water  budget 

•  Plot  head  contours  from  available  well  observations 

•  Sketch  in  flow  lines 

•  Use  simple  radial  well  solutions  (Theis,  Jacob,  etc.)  to  look  at  drawdown 
near  well  fields 

The  desk- top  analysis  should  be  used  to  determine  first  whether  a  modeling 
study  is  necessary  and  then  how  it  should  be  conducted. 

4.  Identify  the  inputs  which  will  have  the  most  effect  on  the  model’s  predictions. 
If  drawdown  is  of  primary  interest  these  will  be  the  aquifer  parameters  and 
the  postulated  incremental  pumpage  and  recharge  rates.  Special  effort  should 
be  devoted  to  the  collection  and  analysis  of  data  needed  to  estimate  these 
parameters. 

5.  Carry  out  a  steady-state  regional  calibration  whenever  possible.  Do  not  rely 
solely  on  localized  pump  test  analyses  when  estimating  spatially-averaged  model 
inputs.  Instead,  use  pump  test  results  to  constrain  a  regional  calibration. 


6,  Make  a  serious  attempt  to  quantify  the  accuracy  of  the  model’s  predictions 
throughout  the  entire  prediction  period.  The  accuracy  evaluation  may  be 
based  on  a  goodness-of-fit  analysis,  a  sensitivity  analysis,  or  both.  Avoid  the 
traditional  “model  validation”  approach. 

7.  When  presenting  model  results,  carefully  consider  the  needs  of  the  reader. 
Are  the  results  presented  in  a  way  which  will  make  decisions  more  informed? 
Are  prediction  errors  acknowledged  or  is  the  accuracy  issue  avoided?  Proper 
presentation  is  crucial  if  the  model’s  results  are  ever  to  be  used. 

The  above  list  is  somewhat  biased  and  probably  not  complete  but  it  does  reflect 
the  major  lessons  learned  in  the  San  Andres-Glorieta  case  study.  If  each  of  these 
guidelines  is  carefully  followed  the  resulting  modeling  study  is  much  more  likely  to 
be  technically  defensible  and  convincing.  The  most  important  factors  influencing 
the  success  of  any  modeling  effort  arc,  of  course,  the  modeler’s  competence  and 
judgment.  There  are  still  no  computerized  substitutes  for  these. 
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